Defect reconstruction is essential in nondestructive testing and structural health monitoring with guided ultrasonic waves. This paper presents an algorithm for reconstructing notches in steel plates, which can be seen as artificial defects representing cracks by comparing measured results with those from a simulation model. The model contains a parameterized notch, and its geometrical parameters are to be reconstructed. While the algorithm is formulated and presented in a general notation, a special case of guided wave propagation is used to investigate one of the simplest possible simulation models that discretizes only the cross section of the steel plate. An efficient simulation model of the plate cross section is obtained by the semianalytical scaled boundary finite element method. The reconstruction algorithm applied is gradientbased, and algorithmic differentiation calculates the gradient. The dedicated experimental setup excites nearly plane wave fronts propagating orthogonal to the notch. A scanning laser Doppler vibrometer records the velocity field at certain points on the plate surface as input to the reconstruction algorithm. Using two plates with notches of different depths, it is demonstrated that accurate geometry reconstruction is possible.
I. INTRODUCTION
Characterizing damage in shelllike structures, such as plates, laminates, or pipes, is an important issue in nondestructive testing and structural health monitoring. Methods based on ultrasonic guided waves are appealing for damage characterization because the guided wave propagation allows the inspection of large areas of the structure from only a few sensor positions. However, determining the damage size using these methods can be challenging, as the various dispersive modes of guided waves do not always permit a straightforward correlation between reflections from damage and damage size.^{1–4} In addition, the wave packet reflected from the defect must be distinguished from the natural wave propagation in the pristine structure for many methods. Therefore, many stateoftheart methods are based on residual signals that compare a measurement of the current, possibly damaged, condition with a baseline measurement of the pristine condition.^{5–7} In order to use these methods, the baseline measurement of the pristine condition must be available, and it must be guaranteed that the pristine condition is free of damage. This requirement can be restrictive, and baselinefree alternatives are desirable.
A possible alternative and baselinefree approach to solve this problem is an algorithm that compares the measured signals with those generated by a simulation model with damage of a specific size. If the simulated signal is close to the measured signal, the damage in the simulation model is a reconstruction of the actual damage.
The main objective of this paper is to validate an algorithm that is an extension of the concept published in Ref. 8, where the algorithm is tested with artificial data from simulation. The algorithm is able to reconstruct a defect defined by certain geometric parameters. The proposed algorithm is derived by solving an inverse problem formulated as an optimization problem. A forward model is iteratively updated until the output matches the measured signals as closely as possible.
The concept of the algorithm has two key aspects for efficient inversion. The first key aspect of the concept is the application of an efficient forward model. The forward model is based on the semianalytic scaled boundary finite element method (SBFEM), first introduced by Song and Wolf,^{9–11} which results in a comparatively small number of degrees of freedom in the final system of equations. The second key aspect of the concept is the gradientbased optimization procedure. The gradient is computed by applying algorithmic differentiation (AD)^{12,13} to the forward model. By computing the gradient using AD, the optimization is fast and accurate.
This paper discusses how the algorithm can be extended to incorporate data from experiments and situations where the excitation has some undetermined properties. Here, the reconstruction of rectangular notches in steel plates using planar ultrasonic guided waves is demonstrated. Investigations of notches as a surrogate for natural cracks in guided wave studies have a long tradition.^{1–4,14,15} While the description of the proposed algorithm can be easily adapted to general problems, this artificial defect is used because it can be manufactured with high precision. The wave propagation and the type of notches investigated in the current work can also be mathematically reduced to a twodimensional (2D) problem by considering only the cross section of the plate. The crosssectional model is well studied in the guided wave literature (e.g., Refs. 1 and 2) and is a simple case with important features of guided wave propagation, such as dispersion and mode conversion. Reducing the mathematical model to the cross section often allows much faster computation, and the problem is more limited in the number of parameters. Consequently, specimens where the waves propagate in only one direction are often used for the first steps of testing an inverse method (e.g., Ref. 16). On the other hand, the experiment must be designed to fit the mathematical model of propagating waves. This paper, therefore, discusses the challenges of validating 2D models through experimentation.
Another favorable condition for the modeling of the experimental setup is a scanning laser Doppler vibrometer (LDV). This complex measurement system is currently required to make a qualitative comparison between the wave shape of the experiment and a simulation model. The advantage is that the model requires no calibration parameters other than the material parameters of the plate. A measurement with a LDV facilitates modeling. In contrast to a measurement with piezoelectric sensors, a LDV measurement can be treated as a simple point evaluation of the velocity, which does not affect the ultrasonic wave itself. Note that the performance of piezoelectric sensors can be highly dependent on the adhesive layer between the sensor and the plate.^{17,18} In Ref. 17, for example, significant changes are observed when the adhesive layer thickness varies in the range from 10 μm to 40 μm. In addition, the performance of the sensorglue combination is temperature dependent (e.g., Ref. 19), although this temperature dependence is only significant for metal plates if there are large temperature fluctuations. Similar concerns are discussed in the qualitative comparison between the experiment and a simulation in Ref. 20. In order to create a simple simulation model that does not require any calibration parameters other than the material parameters of the plate, the excitation sensor and the adhesive layer are not directly modeled, but are replaced by a transfer functionweighted traction. This transfer functionweighted traction is the major extension of the algorithm required for the current experimental validation, because it automatically adjusts the ratio of the guided wave modes excited by the sensor to the measurement. The adjustment ensures independence from the exact parameters of the sensor and the adhesive layer. It should be emphasized that the validation of 2D models is a crucial step towards a general defect reconstruction method, and this example shows that geometric optimization is possible under laboratory conditions.
As an alternative to geometric defect optimization in the context of modelbased reconstruction, defects can be modeled as material changes. This approach is known as full waveform inversion (FWI) and was pioneered for seismic waves.^{21} References 22–24 are examples of FWI for guided waves that include experimental validation. A similar idea can be found in Ref. 25, where wave velocities define the model parameters. Although these results are very promising, the methods require a numerical model that allows for a material change at each point where a defect might exist, making these models computationally expensive. In contrast, the method presented here takes advantage of the fact that large parts of the structure can be calculated semianalytically by assuming a homogeneous material model, resulting in more efficient models.
As an alternative to gradientbased optimization of defect parameters, Ref. 16 presents statistical optimization of surrogate cracks in beamlike structures, while Ref. 26 considers pipes. The presented methods are very versatile and allow the reconstruction of multiple cracks, but they require a large number of simulations to obtain the statistical values. Both references consider onedimensional (1D) models, where the cracks are approximated by special elements.
This paper is organized as follows. Section II summarizes the measurement setup for the two investigated plates. Section III defines the defect reconstruction algorithm and shows which extensions to previous work had to be made to adapt to the experiment. Section IV analyzes the performance of the algorithm and shows the reconstruction results. The conclusion is given in Sec. V. The Appendix summarizes the details of the SBFEM derivation.
II. EXPERIMENTAL SETUP AND DATA ACQUISITION
The reconstruction is studied on two $ 350 \u2009 mm \xd7 350 \u2009 mm \xd7 2 \u2009 mm$ steel plates having a notch near the center. The notches are $ 0.5 \u2009 mm$ wide and $ 300 \u2009 mm$ long, while the notch depths differ with $ 0.4 \u2009 mm$ and $ 0.8 \u2009 mm$, respectively. Figure 1 shows photos of the plates and an overview of the experimental setup, while Fig. 2 presents the schematic layout of the plates. The positioning of the origin of the coordinate system is arbitrary, but for convenience, the origin was chosen to be at the left edge of the notch. This makes deviations of the reconstruction from the nominal notch position trivial to calculate.
As mentioned in the Introduction, the reconstruction algorithm is based on a 2D crosssectional model. However, a crosssectional model simulates only plane wave fronts, which are constant orthogonal to the direction of wave propagation. Therefore, the experimental setup should favor the propagation of plane waves to achieve the best possible agreement between the simulation model and the measurement.
In summary, the data collection process can be divided into the following steps, which will be discussed in more detail below:

The excitation signal is defined on a personal computer (PC).

Based on step 1, the electrical input signal for the sensor is generated and amplified.

The sensor uses the electrical signal to excite the ultrasonic waves inside the plate.

A LDV measures the ultrasonic waves across a scanning area (Fig. 1).

The measurement is postprocessed to create the input for the reconstruction algorithm.
Step 1: Data acquisition begins with the definition of the normalized excitation signal s. Here, the normalized excitation signal is a sinusoidal Gaussian pulse, i.e., $ s ( t ) = sin ( 2 \pi f c t ) \u2009 exp ( \u2212 0.5 ( t \u2212 t c ) 2 / f c \u2212 2 )$ with a center frequency, f_{c}, of $ 500 \u2009 kHz$ and a time shift, t_{c}, of $ 5 f c \u2212 1 = 10 \u2009 \mu s$. The choice of a Gaussian pulse is intended to promote agreement between experiment and simulation, since numerical models for smooth excitations achieve higher accuracy with the same number of degrees of freedom. Figure 3(a) shows the normalized excitation signal in the time domain, while Fig. 3(b) shows the amplitude in the frequency domain with renormalization. In addition, Fig. 3(b) shows the wavelengths of the S0, A0, and S1 modes of the plate. For the corresponding spectrum, only the S0 and A0 modes are excitable. The dispersion relationship between frequency and wavelength is calculated by the SBFEM.^{27,28}
Steps 2 and 3: The normalized excitation signal s is converted into an electrical signal using the arbitrary waveform generator (33521A; Agilent/dataTec GmbH, Reutlingen, Germany) and then amplified by a highvoltage amplifier (HVA400A; Ciprian, Grenoble, France). The amplified signal is sent to a macro fiber composite (MFC) sensor (M8503P2MFC; Smart Material GmbH, Dresden, Germany). The MFC sensor has an active piezoelectric area of $ 3 \u2009 mm \xd7 85 \u2009 mm$, and a total area of $ 7 \u2009 mm \xd7 94.5 \u2009 mm$. The long extension of the sensor in the ydirection compared to the wavelength of the propagating modes (about $ \lambda \u0302 A 0 = 4.6 \u2009 mm$ and $ \lambda \u0302 S 0 = 10.4 \u2009 mm$) produces a nearly planar wavefront and thus to the wave type that also fits the model assumption of a plane strain state. The velocities in the ydirection are decoupled from the velocities in the x and zdirections for plane waves propagating in the xdirection. Due to the choice of the sensor, shear horizontal modes associated with oscillations in the ydirection are negligibly excited, so only velocities in the x and zdirections are considered.
Step 4: A scanning LDV (PSV5003DHV; Polytec GmbH, Waldbronn, Germany) records the resulting ultrasonic waves. Due to the high reflectivity of the steel plate, it was necessary to treat the scanning area with white, nonpermanent paint (ARDROX 9D1B; BASF, Ludwigshafen, Germany) to ensure a defused reflection of the laser (see Fig. 1). The scanning LDV measures the velocity field in all three directions in a predefined grid. For our application, the yvelocity field is neglected. To increase the signaltonoise ratio, the signal is averaged over 1000 times for each measurement point. For this feasibility study, it was decided to test the algorithms on a very high data quality, in most cases 200 replicates are enough for metals (e.g., Ref. 20).
As mentioned above, the model simulates planar waves. Therefore, an evaluation of the planar waves from the experimental data is necessary for a good fit. Figure 4 shows that the excitation already produces a nearly planar wave. Nevertheless, some small changes in the ydirection can be observed. As a filter for waves with a wave vector not parallel to the x axis, the average of the velocities along the y axis is taken. In addition, Fig. 4 shows that the sensor setup excites the S0 and A0 modes simultaneously, which can be distinguished by their wavelength.
Step 5: Since the automatic focusing of the lasers and the alignment of the plate lead to slight deviations from an ideal uniform grid, the x and zvelocity fields are linearly interpolated to a smaller uniform grid. As Fig. 2 shows, the uniform grid has almost the same grid width and is only used to correct the angles and to reduce the number of measurement points. If the exact excitation is known, a measurement at a single xcoordinate is theoretically sufficient to reconstruct the defect parameters^{8} (compare also Ref. 16) However, since in our case the exact excitation is not sufficiently determined, we need several measuring points with different xcoordinates to additionally calculate the transfer function [Eqs. (8) and (9)] of the sensor. The reduced number of measuring points still allows a sufficiently accurate derivation of the transfer function, but also shows that the measuring points can be placed at a significant distance from the defect. The large size of the original scan grid is only required for Fig. 4. In addition, the new uniform grid enables averaging of the data in the ydirection, which, as mentioned above, comes closer to a plane strain signal. This results in a preliminary array $ V \u0303 meas$ of size $ 2 \xd7 11 \xd7 9 \xd7 782$: here, the first dimension distinguishes between the x and zvelocities, the second dimension is for the change in the xcoordinate, the third dimension is for the change in the ycoordinate, and the fourth dimension is for the change in time.
III. THE INVERSE PROBLEM
A. The forward operator
The forward operator mainly calculates the wave propagation in a crosssectional model of a steel plate. The main idea is to use a semianalytical model based on the SBFEM, which leads to a comparatively small number of degrees of freedom in the final system of equations. The SBFEM approximates the solution of the wave equations by dividing the domain into smaller subdomains called superelements. The solution approach for each superelement is derived by a coordinate transformation that uses a finite element approximation of parts of the boundary. The current work uses the (complex) frequency domain to solve the elastic wave equation. For the frequency domain, the finite element boundary approximation contains all degrees of freedom of the approximation, which leads to an order of magnitude reduction in the degrees of freedom compared to a pure finite element approximation of the entire domain. In this work, two types of superelements are considered: on the one hand, the highorder approach for polygonal elements,^{32} and on the other hand, the waveguide superelements.^{33,34}
As mentioned above, the problem is solved in the (complex) frequency domain to obtain a model with very few degrees of freedom. The exponential window method (EWM)^{8,35} is used to compute a finite domain in the frequency domain. The EWM is based on the discrete Laplace transformation, $ DLT \zeta ( s ) = FFT ( s \u2299 \u2009 exp ( \u2212 \zeta t ) )$, and its inverse, $ DLT \zeta \u2212 1 ( s \u0302 ) = FFT \u2212 1 ( s \u0302 ) \u2299 \u2009 exp ( \zeta t )$, where FFT is the fast Fourier transform over the time dimension, s and $ s \u0302$ are arbitrary signal arrays in the time and frequency domains, t is the uniform time vector of the measurement, $ \zeta = 0.5 \omega \Delta $ is the namegiving exponential factor, $ \omega \Delta $ is the resulting angular frequency step due to the sampling frequency of the time vector t, $ exp$ is the entrywise exponential function, and $\u2299$ is the entrywise product along the time dimension. In addition to using these transforms, the wave equation in the EWM is derived for discrete complex angular frequencies: $ \omega \u2113 = ( \u2113 \omega \Delta \u2212 i \zeta ) \u2208 \u2102$. Note that only the real part of the complex angular frequency changes with the index, while the imaginary part remains constant.
For the SBFEM forward model (Fig. 6), the domain is subdivided into superelements. For the current forward model, two different types of superelements are used to calculate a local stiffness matrix, $ S \u0302 b$. The approximation approaches for both types, the waveguide superelement and the polygonal superelement, are briefly summarized in the Appendix. The waveguide superelement is used for the long, undisturbed parts of the crosssection model, as it only requires degrees of freedom in the zdirection and can be of any length (see the Appendix). The polygonal superelement, on the other hand, models all detailed areas of the waveguide. See the detailed views in Fig. 6. It is worth noting that no grid refinement is required for the stress singularities in the notch corners, as these are handled by the semianalytical nature of the approximation approach (see the Appendix).
Note that in the case under investigation, the sensor and the defect are located on different sides of the measurement grid. The derived transfer function can only reproduce the waves propagating in the positive xdirection and not the signals reflected by the defect propagating in the negative xdirection. This observation allows the transfer function and the defect parameters to be correctly determined at the same time. To create a similar situation where the position of the defect can be arbitrarily located, the LDV measurements could be performed on both sides of the sensor. Note that a defect between the sensor and the measurement points could lead to a mathematical ambiguity.
As mentioned in the Introduction, the forward operator in Ref. 8 is mainly extended by the calculation of the transfer function. The number of vectors given by the linear system of Eqs. (7) increases by a factor of 2, but no new square matrices corresponding to all degrees of freedom must be calculated. The 2 × 2 linear system of Eq. (8) and the matrix vector product in Eq. (9) add negligible computational complexity.
A simplified version of the forward operator is shown in Algorithm 1 with the matrices for each the superelement from the Appendix. This version is intended as a guide for a first implementation. However, many steps can be taken to speed up this version. For example, quantities that are independent of the geometric parameters q, such as F, A, $ \omega , \u2009 and \u2009 V \u0302 meas$, can be calculated only once and treated as additional input values. In our implementation, the frequencyindependent matrices, e.g., $ M 0 , \u2009 E 0 , \u2009 E 1 , \u2009 E 2$, are stored for all superelements that are independent of the geometric parameters q. These superelements are to the left of $ x = \u2212 50 \u2009 mm$.
B. The gradient descent and the initial value
The reconstruction algorithm utilizes gradientbased optimization strategies. In contrast to typical derivativefree optimization methods, gradientbased optimization methods have provable convergence properties such as guaranteed local convergence with high convergence speed. However, such algorithms determine only local minima. Since the solution of the parameter estimation problem considered here is given by the global minimizer, for gradientbased methods it is important to provide a good initial guess which is in proximity of the globally optimal point. We achieve this goal by computing the error between simulation and measurement at 100 random points uniformly distributed throughout the parameter space $Q$. The point with the smallest error is selected as an initial guess for the optimization procedure.
As an optimization approach, we selected an iteratively regularized Gauss–Newton method (IRGNM),^{36} as shown in Algorithm 2.
In each iteration n, the IRGNM minimizes a linearization of the original nonlinear leastsquares problem (the error between simulation and measurement) with an additional Thikonov regularization term that vanishes as $ n \u2192 \u221e$. The method is well understood and analyzed and is predominantly used in illposed parameter identification problems. Classical theory also guarantees convergence results in the presence of measurement errors. Different from in prior cases,^{8} when using synthetic data instead of real measurement data, it was not necessary to utilize regularization. Hence, we set $ \alpha n = 0$ for the current work.
This method is highly dependent on accurate and cheaply available gradient information. Derivatives could be computed numerically by applying a finite difference scheme to the forward operator. However, finite difference approximations of the derivative are usually very error prone and may be very costly to compute if the function to be differentiated has many variables. In the past, we have experienced some issues with accuracy in this context and thus chose to use a different approach: As the whole simulation code is available to us, we apply AD.^{12,13} AD as used here is a technique to automatically convert the code of programmed functions with input parameters into code that calculates both the function and the derivative of the function with respect to a given subset of input parameters. The basic principle of the conversion is to decompose the code into its basic commands and insert new lines of code that calculate the derivation with respect to the given parameters of the basic commands using the chain rule. We selected the AD tool adimat (see Ref. 37) in forward mode as it can differentiate simulations written in matlab. Furthermore, adimat can handle some technical programming issues such as data structures, etc., that other AD tools for matlab cannot. Overall, an examination of the computed derivative values showed that the derivatives were within negligible tolerance of expected derivative values. A matlab implementation of an algorithm similar to Algorithm 1 is then translated by adimat into a version that computes $ F ( q )$ and the gradient $ F \u2032 ( q )$ with respect to the defect parameters q.
IV. EVALUATION OF THE PROPOSED ALGORITHM
In this section, we present the performance of our reconstruction algorithm for the two studied plates with different notch depths. We will refer to the plate with a nominal notch depth of $ 0.4 \u2009 mm$ as plate I, while the plate with a nominal notch depth of $ 0.8 \u2009 mm$ will be referred to as plate II. The reconstruction is performed for a single parameter space that includes both notch depths. The position of the notch q_{1} is in the range between −40 mm and $ + 40 \u2009 mm$, which means that the distance is more than 17 A0wavelengths. The notch depth is between $ 0.1 \u2009 mm$ and $ 1.1 \u2009 mm$. The notch width is set to $ 0.5 \u2009 mm$, the nominal width of both plates. The notch width is not of much interest because the notch is only used as a replacement geometry for a crack. Normally, cracks are assumed to be infinitesimally wide, so no width needs to be reconstructed. It is essential that the numerical model is stable for the entire parameter space, since the initial values of the experiments are randomly drawn from this parameter space, as described in Sec. III B.
As a preliminary investigation to the optimization, the objective function is examined. This additional investigation only provides insight into the optimization problem but is not necessary for the actual reconstruction. Figure 7 shows the objective function for both plates over the entire examined parameter space. The colored contour plots are based on a parameter grid of size 161 × 21. Both plates show a complex shape for the objective function with several local minima. However, the global minimum is located near the nominal parameters $ q * = ( 0 \u2009 mm , \u2009 0.4 \u2009 mm ) T$ and $ q * = ( 0 \u2009 mm , \u2009 0.8 \u2009 mm ) T$, respectively. There are several reasons why the global minimum of the objective function is not exactly at the nominal parameters. First, the measurement naturally contains a certain amount of noise. Also, the assumptions for the numerical model may not be accurate enough. It should be emphasized that this is only a 2D crosssectional model. In addition, the plates were machined with a certain tolerance so that the actual values differ from the nominal parameters. To illustrate the machiningrelated variance, the statistical values of several gauging points at which the notch is assessed with a measuring microscope (MS4BM01; Uhl Technische Mikroskopie, Aßlar, Germany) are listed in Table I. The distances were measured at 14 equally distributed points between $ y = \u2212 130 \u2009 mm$ and $ y = + 130 \u2009 mm$ along the notches (see Fig. 2). From the minimum and maximum depth values, it is evident that the bottoms of the notches are not completely flat, but there is no clear trend from one side of the notch to the other. On average, the notches are deeper than the nominal values. For the sake of simplicity, and since q_{1} varies over the length of the notch, the nominal parameter is usually referred to in the following. Furthermore, the MFC sensors were manually aligned to the predefined position and glued on by hand. The measuring points were aligned visually using the camera system of the LDV, so a more precise determination than the nominal position of $ q 1$ cannot be given. Nevertheless, Fig. 7 shows that the reconstruction parameters in the area with a low objective function value are close to the nominal parameters.
Plate .  Width (mm) .  Depth (mm) .  

Minimum .  Mean .  Maximum .  Minimum .  Mean .  Maximum .  
Plate I  0.498  0.504  0.507  0.408  0.427  0.450 
Plate II  0.496  0.500  0.505  0.763  0.809  0.850 
Plate .  Width (mm) .  Depth (mm) .  

Minimum .  Mean .  Maximum .  Minimum .  Mean .  Maximum .  
Plate I  0.498  0.504  0.507  0.408  0.427  0.450 
Plate II  0.496  0.500  0.505  0.763  0.809  0.850 
Gauged at 14 positions between $ y = \u2212 130 \u2009 mm$ and $ y = + 130 \u2009 mm$ (see Fig. 2).
The final reconstruction with the corresponding parameters $ q min$ (Table II) will be examined in the following. The initial parameter vectors are $ q 0 = ( 0.24 \u2009 mm , \u2009 0.34 \u2009 mm )$ and $ q 0 = ( \u2212 0.07 \u2009 mm , \u2009 0.59 \u2009 mm )$ for plates I and II, respectively. The final reconstruction is achieved after 17 and 13 iterations, where ϵ in Algorithm 2 is set to 1 × 1^{−5}. Figure 8 shows the spectrum of the calculated transfer functions $ h \u0302 i$ [Eq. (8)]. It can be seen that both tractions make a significant contribution to the signal. Moreover, the maximum value of the transfer functions is close to the center frequency of $ 500 \u2009 kHz$ of the input signal $ s \u0302$. However, there are also qualitative differences between the two cases. While Fig. 8(a) shows an excellent agreement in the overall shape of the spectra between the transfer functions $ h \u0302 i$ and the input signal $ s \u0302$, a clear difference in shape can be observed in Fig. 8(b). The transfer functions represent not only the different MFC sensors, but also the adhesive bond to the plate. The fact that the adhesive bond is handmade could explain the differences, despite the same design of the sensors and exactly the same signal generator and amplifier. However, the difference shows that calculating the transfer functions makes the optimization more robust to variations in sensor performance.
Plate .  Parameters .  Measured mean depth (mm) .  

Nominal .  Reconstruction .  
Plate I  $ q * = ( 0 \u2009 mm , \u2009 0.4 \u2009 mm ) T$  $ q min = ( \u2212 0.09 \u2009 mm , \u2009 + 0.44 \u2009 mm ) T$  $ q 2 meas = 0.427 \u2009 mm$ 
Plate II  $ q * = ( 0 \u2009 mm , \u2009 0.8 \u2009 mm ) T$  $ q min = ( + 0.03 \u2009 mm , \u2009 + 0.77 \u2009 mm ) T$  $ q 2 meas = 0.809 \u2009 mm$ 
Plate .  Parameters .  Measured mean depth (mm) .  

Nominal .  Reconstruction .  
Plate I  $ q * = ( 0 \u2009 mm , \u2009 0.4 \u2009 mm ) T$  $ q min = ( \u2212 0.09 \u2009 mm , \u2009 + 0.44 \u2009 mm ) T$  $ q 2 meas = 0.427 \u2009 mm$ 
Plate II  $ q * = ( 0 \u2009 mm , \u2009 0.8 \u2009 mm ) T$  $ q min = ( + 0.03 \u2009 mm , \u2009 + 0.77 \u2009 mm ) T$  $ q 2 meas = 0.809 \u2009 mm$ 
Input: geometric parameter q, time vector t, the measurement array $ V meas$, regularization parameter β . 

Block of all frequencyindependent quantities: 
1: compute ζ from t 
2: $ V \u0302 meas = DLT \zeta ( V meas )$ 
3: allocate $ V \u0302 sim$, A, F 
4: compute A 
5: for all finite elements e do 
6: compute $ F e$ 
7: add $ F e$ to F 
8: for all superelements do 
9: compute and store $ M 0 , \u2009 E 0 , \u2009 E 1 , \u2009 E 2$ depending on the type of the superelement 
10: if the current superelement is a polygonal superelement then 
11: for m = 1,…,M do 
12: compute and store $ S 0 ( m ) , \u2009 S 1 ( m ) , \u2009 X ( m + 1 )$ 
Block of all frequencydependent quantities: 
13: determine a relevant frequency range $\omega $ from $ V \u0302 meas$ 
14: for all ω in $\omega $ do 
15: allocate $ S \u0302$ 
16: for all superelements do 
17: compute $ S \u0302 b$ depending on the type of the superelement 
18: add $ S \u0302 b$ to $ S \u0302$ 
19: $ U \u0302 = S \u0302 \ F$ 
20: $ V \u0302 sim = ( i \omega ) A U \u0302$ 
21: $ v \u0302 meas = vec ( V \u0302 meas  \omega )$ 
22: $ h \u0302 = ( V \u0302 sim H V \u0302 sim + \beta I ) \ ( V \u0302 sim H v \u0302 meas )$ 
23: $ v \u0302 sim = V \u0302 sim h \u0302$ 
24: sort $ v \u0302 sim$ into $ V \u0302 sim$ 
Transform back to time domain and final steps: 
25: $ V sim = DLT \zeta \u2212 1 ( V \u0302 sim )$ 
26: $ y sim = vec ( env t ( V sim ) )$ 
Input: geometric parameter q, time vector t, the measurement array $ V meas$, regularization parameter β . 

Block of all frequencyindependent quantities: 
1: compute ζ from t 
2: $ V \u0302 meas = DLT \zeta ( V meas )$ 
3: allocate $ V \u0302 sim$, A, F 
4: compute A 
5: for all finite elements e do 
6: compute $ F e$ 
7: add $ F e$ to F 
8: for all superelements do 
9: compute and store $ M 0 , \u2009 E 0 , \u2009 E 1 , \u2009 E 2$ depending on the type of the superelement 
10: if the current superelement is a polygonal superelement then 
11: for m = 1,…,M do 
12: compute and store $ S 0 ( m ) , \u2009 S 1 ( m ) , \u2009 X ( m + 1 )$ 
Block of all frequencydependent quantities: 
13: determine a relevant frequency range $\omega $ from $ V \u0302 meas$ 
14: for all ω in $\omega $ do 
15: allocate $ S \u0302$ 
16: for all superelements do 
17: compute $ S \u0302 b$ depending on the type of the superelement 
18: add $ S \u0302 b$ to $ S \u0302$ 
19: $ U \u0302 = S \u0302 \ F$ 
20: $ V \u0302 sim = ( i \omega ) A U \u0302$ 
21: $ v \u0302 meas = vec ( V \u0302 meas  \omega )$ 
22: $ h \u0302 = ( V \u0302 sim H V \u0302 sim + \beta I ) \ ( V \u0302 sim H v \u0302 meas )$ 
23: $ v \u0302 sim = V \u0302 sim h \u0302$ 
24: sort $ v \u0302 sim$ into $ V \u0302 sim$ 
Transform back to time domain and final steps: 
25: $ V sim = DLT \zeta \u2212 1 ( V \u0302 sim )$ 
26: $ y sim = vec ( env t ( V sim ) )$ 
Input: initial value $ q 0$, regularization parameter sequence $ ( \alpha n ) n \u2208 \mathbb{N} \u2192 0$, maximal number of iterations N . 

1: for $ n = 0 , \u2026 , N$ do 
2: $ q n + 1 = q n + ( F \u2032 ( q n ) H F \u2032 ( q n ) + \alpha n I ) \u2212 1 ( F \u2032 ( q n ) H ( y meas \u2212 F ( q n ) ) + \alpha n ( q 0 \u2212 q n ) )$ 
3: if $ \u2225 q n + 1 \u2212 q n \u2225 < \u03f5$, then 
4: $ q min = q n + 1$ 
5: STOP 
Input: initial value $ q 0$, regularization parameter sequence $ ( \alpha n ) n \u2208 \mathbb{N} \u2192 0$, maximal number of iterations N . 

1: for $ n = 0 , \u2026 , N$ do 
2: $ q n + 1 = q n + ( F \u2032 ( q n ) H F \u2032 ( q n ) + \alpha n I ) \u2212 1 ( F \u2032 ( q n ) H ( y meas \u2212 F ( q n ) ) + \alpha n ( q 0 \u2212 q n ) )$ 
3: if $ \u2225 q n + 1 \u2212 q n \u2225 < \u03f5$, then 
4: $ q min = q n + 1$ 
5: STOP 
The optimization procedure should achieve the best possible agreement between simulation and measurement at the measurement points. To give an impression of the accuracy of the model, Fig. 9 shows a direct comparison of the signals. Although the inverse method only uses the envelopes to compute the objective function, the phase of the signal from the simulation also matches the measured signal. The matching phase can be explained by the calculation of the transfer function, where the phase information is adjusted.
For the evaluation of the reconstruction, we refer to Fig. 10. In this figure, the SBFEM mesh of the reconstruction is plotted against the nominal geometry outlined by a dashed line. The figures show only a small section around the notch. It should be noted that the actual parameter space is much larger. The figure shows a typically numbered x axis at the bottom, while an x axis measured with the A0wavelength of the center frequency f_{c} is shown at the top. This additional axis should make it possible to estimate the error for other plate dimensions, since the wavelength is the characteristic length of the problem. The errors for the xposition $  q 1 * \u2212 q 1 min $ of the notch are $ + 0.09 \u2009 mm$ and $ + 0.03 \u2009 mm$ for plates I and II, respectively. In relation, the nominal errors for the notch depth $  q 2 * \u2212 q 2 min $ are $ + 0.04 \u2009 mm$ and $ + 0.03 \u2009 mm$, while the errors to average measured notch depth $  q 2 meas \u2212 q 2 min $ are $ 0.01 \u2009 mm$ and $ 0.04 \u2009 mm$.
V. CONCLUSION
An algorithm has been presented that enables defect reconstruction using guided ultrasonic waves by solving an inverse optimization problem. The output of a simulation model, which contains a parametric representation of the defect, is intended to approximate the experimental data as closely as possible. For this purpose, a gradientbased method is used, where the computation of the derivatives according to the defect parameters is based on AD. A second aspect of the presented work is the highly efficient ABFEM, which requires very small, final linear systems of equations to approximate the linear waveguide equation.
This article presents the successful reconstruction of two notches in steel plates based on data from experiments. Although the reconstruction algorithm is formulated in a very general way, the initial focus was on the simplest case of wave propagation, where the ultrasonic wave can be assumed to be planar, and the problem can be reduced to a 2D crosssectional model. An error between nominal geometry and reconstruction on the order of 0.06 mm was achieved. Thus, it could be shown that, first, a geometric reconstruction of defects is possible and, second, a crosssectional model is capable of describing the wave field with sufficient accuracy for a reconstruction, given a suitable experimental setup.
The first important finding of this work is that the envelopes of the recorded data provide sufficient information while leading to an objective function that allows optimization with a gradient method. It should be noted that the current sensor setup leads to mixed mode excitation, and both modes are converted into each other at the defect. In addition, the modes are subject to dispersion. Nevertheless, the notch parameters are correctly determined by the approach because these effects are taken into account within the model.
The second finding is that another linear fit must be introduced for the sensor's force distribution. This allows reconstruction even for an ultrasound source that is not fully characterized, where the ultrasound source is determined by both the sensor and the adhesive layer.^{17,18} However, realworld applications would benefit from a sensor model that includes the adhesive layer and allows quantitative comparison between simulation and measurement with a small error—eliminating the need for LDV measurement.
Future work will extend the reconstruction algorithm to threedimensional (3D) models of wave propagation. Reconstructing a crack in three dimensions introduces the y position of the crack, the angle of the crack in the xy plane, and the crack length as important parameters. This leads to new challenges. For example, the space of possible initial values is at least fivedimensional (5D), so a random draw is unlikely to come close to the true crack parameters. In addition, if there is only one source of ultrasonic waves, the crack may be orthogonal to the wavefront, resulting in very small reflections. This makes a single source technique unlikely to detect all cracks.
ACKNOWLEDGMENTS
The authors gratefully acknowledge the German Research Foundation for funding (DFG Project No. 428590437). In addition, the authors would also like to thank Yevgeniya Lugovtsova and Tobias Homann for their support with the LDV measurement and Andreas Knöppchen for the data obtained with the measuring microscope.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: THE TYPES OF SUPERELEMENTS
This appendix briefly summarizes the derivation of the two types of superelements used and introduces the notation of the matrices for Algorithm 1.