The U.S. Air Force seeks to improve lifecycle management of composite structures. Nondestructive characterization of damage is a key input to this framework. One approach to characterization is model-based inversion of ultrasound inspection data; however, the computational expense of simulating the response from damage represents a major hurdle for practicality. A surrogate forward model with greater computational efficiency and sufficient accuracy is, therefore, critical to enable damage characterization via model-based inversion. In this work, a surrogate model based on Gaussian process regression (GPR) is developed on the chirplet decomposition of the simulated quasi-shear scatter from delamination-like features that form a shadowed region within a representative composite layup. The surrogate model is called in the solution of the inverse problem for the position of the hidden delamination, which is achieved with <0.5% error in <20 min on a workstation computer for two unique test cases. These results demonstrate that solving the inverse problem from the ultrasonic response is tractable for composite impact damage with hidden delaminations.

## Introduction

The increasing use of polymer–matrix composites within Air Force structures has fueled a growing demand for improved lifecycle management of these materials [1]. The Air Force Aircraft Structural Integrity Program guidelines for metallics provide the template for this improvement, where validated damage progression models enable tolerance to assumed defects for an acceptable measure of risk [2]. Composites damage progression models are sensitive to the geometric complexity of typical damage modes. Consequently, nondestructive evaluation (NDE) techniques are needed to volumetrically represent damage for improved lifecycle management.

Currently viable composite NDE methods include X-ray-computed tomography (XCT) and ultrasound. XCT scans volumetrically reconstruct internal damage, making it an obvious path for feature characterization. Nevertheless, XCT cannot be used in situ as a single-sided inspection and generally requires additional safety precautions. By contrast, normal incidence longitudinal wave ultrasound can detect the presence of delaminations and matrix cracking in a single-sided inspection (Fig. 1(a)). However, it is incapable of insonifying the hidden damaged region under the top delamination layers (Fig. 1(b)). A novel angle beam shear wave method addresses this limitation by using direct and back-wall reflected quasi-shear waves to inspect for hidden ply delaminations and matrix cracks [3,4]. This approach builds on prior work on polar backscatter ultrasonic techniques [5–8]. In this way, additional information from hidden damage features can be captured and processed to characterize the damage producing the observed response. It should be noted that phased array transducers (either in a normal incidence or angle beam shear configuration) may provide additional information; however, for the purposes of this work, the simplest sensing approach (single element, pulse-echo inspection) was selected based on the complexity of the material and internal damage.

One challenge with this angle beam ultrasonic inspection technique is resolving the multiple signals scattered from delamination edges at different depths. Manual interpretation of the paths for each signal is complicated, as the scattered signals interact with the composite panel surfaces and neighboring delaminations. Model-based inversion is proposed to address this challenge. This approach uses nonlinear least-squares based optimization to estimate the parameters of a physics-based forward model by minimizing an objective function defined as the squared residual between the model output and the observed response. Prior work has considered flaw reconstructions with ultrasonic NDE in anisotropic materials [9–13]. Other efforts have examined the use of multipath signals for characterizing the location of scatterers [14–17]. For example, Hutt and Simonetti used a strongly scattering planar interface as a mirror to look behind the object to achieve complete reconstruction [17]. Some studies have investigated the identification of multiply diffracted ultrasonic waves from point scatterers, although many of them have focused on medical ultrasound applications [18,19]. Most approaches outlined within the literature have leveraged approximate models like ray theory for improving reconstruction. The proposed solution builds on prior work using numerical solutions of parameterized models to invert the damage state [20–22]. The specific problem of inverting the location of multiple impact damage discontinuities in a thin composite panel has not been addressed to date. In this work, a simplified case study was investigated to demonstrate the model-based inversion approach.

## Materials and Methods

The geometry of impact damage within a composite is dependent on the size, shape, and velocity of the impactor; the interaction of the deforming impactor with the shape and orientation of the structure; and the current mechanical state of the material as a function of as-manufactured properties, internal defects, thermo-mechanical degradation, and environment. These variables result in asymmetric, highly irregular damage contours. Bench-top impact tests on laboratory coupons can control some of these variables, albeit with little reduction in the geometric complexity of the damage.

A simple representative problem, estimation of the location of a hidden delamination beneath a top delamination, was considered to reduce geometric complexity to a level suitable for a feasibility study. Two cases are presented in Fig. 2: in case 1, the lower (hidden) delamination was fixed in *z* with the *x*-position iterated through nine values; in case 2, the lower delamination was fixed in *x* with the *z*-position iterated through eight values.

The parameters for the study follow existing specimens and laboratory testing conditions [23]. A 3.2 mm thick, 24 ply IM7/977-3 composite panel of [−45 deg/90 deg/45 deg/0 deg]_{3 s} layup with simplified delamination-like features was considered. The simulated ultrasonic transducer was 6.35 mm in diameter, with a center frequency of 5 MHz and a focal length of 19.05 mm. The angle of incidence was set to 24 deg, near the first critical angle of the homogenized composite (∼27 deg) such that the incident beam was refracted into quasi-shear modes at the front-wall. The simulated delaminations were each 12.7 mm long. Their locations did not directly coincide with ply interfaces but did not exceed a minimum of a single ply separation in the *z*-direction. These cases provide an elegant demonstration of the proposed method of solving the inverse problem while leaving room for follow-on studies addressing the challenges of a much larger set of pertinent parameters.

## Calculation

### Forward Model and Need for Surrogate Representation.

The efficiency of model-based inversion is proportional to the well-posedness of the problem and inversely proportional to the cost of each forward model iterate. Two computational tools—CIVA FIDEL 2D and PZFlex—are tailored for ultrasound simulation of composite materials [24–26]. CIVA FIDEL 2D couples semi-analytical beam calculations with a finite difference time domain (FDTD) solver acting on a user-specified region to model two-dimensional (2D) wave propagation for simple geometries and boundary conditions [27,28]. Comparatively, PZFlex uses finite element analysis with an optimized explicit solver to simulate 2D or three-dimensional wave propagation on any geometry with complex boundary conditions. Both tools can simulate a spatially dense B-scan response from thin delamination-like features in a representative composite material within several hours; however, typical inverse problems often require hundreds of model calls. For a 2D ultrasonic inspection scenario, solution times would range from days to weeks. When the physics-based model is computationally expensive, a surrogate will often be developed that sacrifices some accuracy for a decrease in computational expense [29]. The need for an efficient surrogate forward model is evident for the ultrasound inverse problem; thus, CIVA was chosen for surrogate model development based on its greater computational efficiency, the dimensionality of the scenario, and previous validation work [4].

Previous work demonstrated the feasibility of an ultrasound surrogate model for the case of a single delamination placed in four unique *z*-positions [30]. In this work, a similar model is developed on the more challenging case of a hidden delamination. Model development began with selection of the representative problem, which was modeled in CIVA FIDEL 2D for a sparse set of design points containing the inverse solution. Chirplet decomposition was applied to each simulated B-scan; then, the surrogate model was constructed by interpolating on the chirplet parameter space via Gaussian process regression (GPR). The resulting fits formed the surrogate model used to solve the inverse problem for the ultrasonic response from a hidden delamination with an unknown *x*- or *z*-position.

### B-Scan Simulation.

Several assumptions were made for the simulation. First, it was assumed that the composite was flat and parallel, which is valid for a simulation-based study but likely invalid for real composite panels. Second, it was assumed that delaminations could be modeled as discontinuities with zero thickness parallel to the plies. While acceptable for a simulation-based study, real delamination morphologies should be incorporated to validate this assumption. Finally, it was assumed that water is present on the backside of the composite. This assumption would hold true for a real panel scanned via immersion ultrasound; however, an air back-wall is more likely for real inspection scenarios.

The model consisted of a water path to simulate an immersion-based inspection, the plies forming the composite panel, and delamination features. Because of the FIDEL package, the plies were modeled individually (the composite was not homogenized). For these boundary conditions, a 2D FDTD region was set to surround the near edges of the delaminations and the composite panel boundaries. The FDTD solution region moves with the scanned probe, so the lateral extent was set to a large size (20 mm) in order to include both direct and full skip interactions with the far-wall. Perfectly matched layers were also applied to the lateral boundaries to help control fictitious signals from the model domain edges. A mesh resolution was fixed to 1/20th of the center frequency wavelength. Additional details of the formulation can be found in published literature [27,28]. B-scan images were generated by scanning the transducer over a 6 mm path across the top of the composite in 0.1 mm increments with a constant water path of 17 mm. Simulations were performed on a HP computer with two 12-core hyperthreaded processers (48 compute cores) at 2.5 GHz clock speed and 98 GB of RAM. Each B-scan required ∼360 min of simulation time. An example simulated B-scan for a single delamination located at 1.0 mm in depth from the top surface is shown in Fig. 3(a). In this case, the largest signal is a direct diffraction from the delamination edge. Signals later in time correspond to the multiple paths for diffracted responses interacting with the top or bottom surfaces of the composite. Half-skips involve a single reflection off of the back-wall surface, while full-skips involve two reflections off of the back-wall. These signals shift in time and space according to the position of the delamination edge.

### Chirplet Decomposition.

where *β* is the amplitude of the Gaussian envelope, *α*_{1} is the bandwidth of the Gaussian envelope, *τ* is the time of arrival of the response, *f _{c}* is the center frequency,

*ϕ*is the phase angle, and

*α*

_{2}is the chirp rate. The chirp rate,

*α*

_{2}, captures the variation of frequency with respect to time and represents a substantial improvement over a wavelet-based approach. The chirplet has been shown to model individual ultrasonic reflections as a sparse, energy preserving representation of the data without any a priori information on the expected response [31]. The chirplet parameters can be estimated from the analytic signal and further refined using a Gauss–Newton algorithm [32]. In the chirplet decomposition algorithm suggested in Refs. [31–33], chirplets are fit successively to individual A-scans until a reconstruction error tolerance is reached. The method used here differs slightly; rather than faithfully reconstructing each A-scan, chirplets only reproduced reflections in the B-scan that represented responses from a delamination. This was achieved by screening out reflections from the front-wall and back-wall of the composite—the remaining reflections were visually segmented and fitted. The surrogate model predicts how these individual reflections change with the position of the hidden delamination. Chirplet decomposition reduces the number of parameters needed to represent the B-scan to $6\u2211j=1Nnj$, where

*N*is the total number of distinct reflections, and $nj$ is the number of transducer positions containing a signal from the

*j*th reflection. As an example, Fig. 3 depicts a B-scan containing ∼80k samples, while the chirplet reconstruction is defined by ∼360 chirplet parameters [34].

### Gaussian Process Regression.

*f*at a test point

*u*is described by

The assumed mean function, *μ*, describes the general trend of the data. Both the covariance vector *c _{u}* between the test point(s)

*u*and the sample points

*s*and the covariance matrix

*C*between the sample points

*s*depend on the assumed covariance function. In the previous formulation, the known responses

*y*recorded at the sample points

*s*are assumed to be the sum of the unknown latent function values

*f*(

*s*) and independent, identically distributed Gaussian noise with variance

*σ*

_{n}^{2}. For this application, the inputs

*s*are given by the transducer location and

*x*- or

*z*-position of the delamination, while the outputs

*y*are the value of the chirplet parameter

*β*,

*α*

_{1},

*τ*,

*f*,

_{c}*ϕ*, or

*α*

_{2}. For the case of simulation data, the variance

*σ*

_{n}^{2}is assumed to be zero. The outputs

*f*are the estimates of the chirplet parameters for the transducer locations at intermediate

*x*- or

*z*-positions of the hidden delamination.

*u*are the test point(s),

*c*is a constant,

*a*are polynomial coefficients, and

_{i}*n*is the order of the term. In Eq. (4),

*s*are the sample points and ǁ ǁ

_{2}denotes the Euclidean distance. The parameter

*l*describes how rapidly the correlation between samples decreases and is a hyperparameter of the GPR model optimized via leave-one-out cross validation on the input data; thus, the value of

*l*varies between parameters and reflections [38]. Selection of the mean and covariance functions is often guided by physics-based expectations of the response; however, mathematically rigorous model selection techniques such as the Akaike information criterion (AIC, Eq. (5)) can be applied

The AIC rewards model accuracy (the first term, where *N* is the number of samples) while penalizing complexity (the second term, where *d* is the number of model parameters), with better models producing more highly negative values [39]. For this work, the response from each parameter was fit with every combination of the selected mean and covariance functions, with the combination producing the lowest AIC selected as the optimal set for that parameter. The surrogate model results from the best performing mean and covariance function combination with an optimized hyperparameter *l*. An input vector *u* is developed for each intermediate *x*- or *z*-position of the hidden delamination and applied to the GPR model. The output of the model is the estimate of the selected parameter for a given reflection.

## Results and Discussion

### Case 1: *x*-Varying, *z*-Fixed.

For this demonstration, the *z*-position of the lower delamination was fixed at 2 mm from the upper delamination, while *x* ranged from 0.8 mm to 2.4 mm in nine steps from the leftmost edge of the upper delamination. This range was selected due to the evident variations in reflection position and intensity at each delamination position. Changes between B-scans diminish drastically beyond *x* = 2.4 mm, resulting in regions where unique solutions to the inverse problem hypothetically do not exist. The temporospatial translations of nine reflections were tracked across each B-scan. These reflections were selected in order of decreasing amplitude, with the bounds of each reflection governed by the envelope of the chirplet in the time dimension and the threshold for fitting a chirplet to an observed signal in the spatial direction. The B-scans used for building the surrogate model are presented in Fig. 4. Tracked reflections generally maintain their original shape as the hidden delamination varies from *x* = 0.8 mm to *x* = 1.4 mm, although the strength of the response from later reflections begins to diminish. From *x* = 1.4 to *x* = 2.0 more reflections appear, while those later in time continue to decrease in amplitude. Finally, later reflections diminish to the point where their amplitude is below the threshold for chirplet fitting between *x* = 2.0 and *x* = 2.4.

Sample data were drawn from Gaussian chirplet parameters for a given delamination *x*-position. Three mean basis functions (constant, linear, and quadratic) (Eq. (3)) and two covariance basis functions (Eq. (4)) were competed for each parameter fit using the AIC. In general, the Gaussian covariance function provided the most accurate representation of the observed data. The optimal mean function varied from parameter to parameter, although less sensitivity to the mean function was observed. To demonstrate the surrogate model for case 1, selected Gaussian chirplet parameters and the GPR model fits from the nine reflections for *x* = 1.3 mm are presented in Fig. 5. Several of the parameters vary smoothly across a given reflection (e.g., amplitude and time of arrival, Figs. 5(a) and 5(c), respectively), although few parameters are uniformly smooth across all reflections (e.g., reflection eight phase and reflection three chirp rate, Figs. 5(e) and 5(f), respectively). As the quality of the chirplet representation is most sensitive to the amplitude and time of arrival of the signal, the smoothness of the response from these parameters is highly encouraging for accurate representation of an intermediate B-scan location.

Fitted parameters were used to reconstruct the B-scan for test positions of *x* = 1.2 mm and *x* = 2.0 mm as presented in Figs. 6(a)–6(c) and 6(d)–6(f), respectively. The quality of the surrogate was judged by the residual between the simulated B-scan for that delamination location and the estimated B-scan. Construction of the B-scan from the surrogate model took 52 ms. As the nine tracked reflections were selected based in order of decreasing amplitude, the residual for both cases is dominated by low-amplitude reflections; however, some response from tracked reflections (most notably the residual arising from the first, largest reflection) was not fully modeled by the chirplet fits.

The true performance of the surrogate model is measured by its ability to rapidly and accurately invert a B-scan with a delamination in an unknown position. This was tested by first simulating a CIVA B-scan for an intermediate point in the surrogate model (lower delamination positioned at *x* = 1.3 mm). Next, an optimization algorithm was applied to minimize the residual sum-of-squares between the supplied B-scan and that generated by the surrogate model. As reflections may appear or disappear as the *x*-position of the hidden delamination is changed, an algorithm that handles discontinuities in the search space was required. The differential evolution (DE) algorithm was, thus, selected [40].

The CIVA B-scan, initial guesses and final solution of the DE algorithm, surrogate model evaluated at the final solution, and the residual in decibels are presented in Figs. 7(a)–7(d), respectively.

Twenty iterations of the DE algorithm were evaluated using 50 initial points covering the possible solution space, resulting in a final solution of *x* = 1.297 mm. This result has an error of 0.22% with a runtime of ∼18 min on a standard workstation. Visually, the B-scan reconstructed for *x* = 1.297 mm appears nearly indistinguishable from the CIVA B-scan at *x* = 1.3 mm, further demonstrating the quality of the inversion.

### Case 2: *z*-Varying, *x*-Fixed.

For case 2, the *x*-position of the lower delamination was fixed at 1.2 mm from the leftmost edge of the upper delamination, while *z* ranged from 1.1 mm to 2.5 mm from the upper delamination in eight steps. As before, the range of hidden delamination *z*-positions was chosen to ensure distinct variations in the intensity and position of the tracked reflections. The temporospatial translations of nine reflections were tracked across each B-scan. Reflections were selected in order of decreasing amplitude, with the bounds of each reflection governed by the envelope of the chirplet in the time dimension and the threshold for fitting a chirplet to an observed signal in the spatial direction. The B-scans developed for the surrogate model are presented in Fig. 8. As the hidden delamination moves from *z* = 1.1 mm to *z* = 1.5 mm, the reflections around 27 *μ*s increase in size and complexity, although the general shape of the response remains mostly constant. Complex reflections beyond 25 *μ*s appear and shift with the *z*-position of the hidden delamination from *z* = 1.5 mm to *z* = 2.1 mm.

Fitted parameters were used to reconstruct the B-scan for test positions of *z* = 1.3 mm and *z* = 2.1 mm as presented in Figs. 9(a)–9(c) and 9(d)–9(f), respectively. The quality of the surrogate was judged by the residual between the simulated B-scan for that delamination location and the estimated B-scan.

Construction of the B-scan from the surrogate model took 61 ms. As in case 1, the residual is dominated by low amplitude reflections or incomplete fitting of the primary reflection.

To test this case, a CIVA B-scan was generated for the lower delamination positioned at *z* = 2.2 mm. The CIVA B-scan, initial guesses and final solution of the DE algorithm, surrogate model evaluated at the inverse solution, and residual in decibels are presented in Figs. 10(a)–10(d), respectively. The DE algorithm was again used with 20 iterations on 50 initial points covering the possible solution space. An inverse solution of *z* = 2.191 was found. This result has error of 0.41% with a runtime of ∼14 min on a standard workstation. As observed in case 1, the B-scan reconstructed for *z* = 2.191 mm appears nearly indistinguishable from the CIVA B-scan at *z* = 2.2 mm, further demonstrating the quality of the inversion.

## Conclusions

The ultrasonic response from a composite material with embedded delamination-like features was simulated for two test cases. Gaussian chirplets were fit to selected reflections within the simulated B-scans. Chirplet parameters and the position of the shifting coordinate (*x* or *z*) were used to develop GPR surrogate models. These models were used to estimate the unknown position of the hidden delamination for each test case. The unknown positions were rapidly estimated to within <0.5% of the actual value, demonstrating that an efficient surrogate model can be developed on ultrasound simulations and that inversion of the ultrasonic response from delamination-like features is tractable for composite impact damage with hidden delaminations.

Practical application begins with an initial model of the material instantiated with the general shape of the delamination field based on a priori knowledge of relevant parameters (impact energy, layup, etc.). B-scans for permutations of this field would be generated, decomposed by the chirplet transform, and parameterized with Kriging to produce a surrogate model. The residual between the parameterized model and the B-scan of the actual damage would be minimized, resulting in damage features matching those within the real composite structure.

Several challenges drive future work. First, chirplet decomposition relies on accurately mapping distinct reflections within the B-scan. While some degree of separation between responses is imposed by layup and ply thickness, the proximity of damage features often results in merged reflections. The current mapping approach requires user interpretation of the total response, which only compounds the challenge of selecting reflections to fit within the surrogate model. Algorithmic identification of individual reflections is, thus, crucial to enabling the inverse solution for complex damage morphologies. Second, the inverse problem was demonstrated for design points where reflections from the hidden delamination remain within the B-scan. These indications vanish as the hidden feature becomes further occluded by the upper delamination. Successive positions of the hidden delamination will eventually result in the same observed response, resulting in nonunique solutions to the inverse problem. Better methods of insonifying the hidden region is thereby critical to extending the range of tractable design cases. Third, experimental noise and variability from factors such as probe misalignment, water path, and transducer representation (e.g., dimensions/properties of the modeled piezoelectric crystal, backing material, or lens) are anticipated to impact the quality of the inverse solution. Parameter sensitivity studies will be employed to understand the impact of these sources of uncertainty. Finally, the techniques developed within this work could be extended to other relevant composite damage features, including fiber breakage and matrix cracking. Each feature type dominates failure for specific composite configurations, rendering them equally important in predicting the remaining loading cycles survivable by the damaged component.

## Acknowledgment

The authors would like to acknowledge the support of Dr. Amanda Criner for guidance and technical support.

## Funding Data

Air Force Research Laboratory Materials and Manufacturing Directorate (Contract Nos. FA8650-14-D-5224 and FA8650-15-D-5231).