## Abstract

The accurate quantification of wall loss caused by corrosion is critical to the reliable life estimation of pipes and pressure vessels. Traditional thickness gauging by scanning a probe is slow and requires access to all points on the surface; this is impractical in many cases as corrosion often occurs where access is restricted, such as beneath supports where water collects. Guided wave tomography presents a solution to this; by transmitting guided waves through the region of interest and exploiting their dispersive nature, it is possible to build up a map of thickness. While the best results have been seen when using the fundamental modes A0 and S0 at low frequency, the complex scattering of the waves causes errors within the reconstruction. It is demonstrated that these lead to an underestimate in wall loss for A0 but an overestimate for S0. Further analysis showed that this error was related to density variation, which was proportional to thickness. It was demonstrated how this could be corrected for in the reconstructions, in many cases resulting in the near-elimination of the error across a range of defects, and greatly improving the accuracy of life estimates from guided wave tomography.

## 1. Introduction

Corrosion presents a major challenge for the petrochemical industry, reducing the wall thickness of pipes and pressure vessels and increasing the chance of rupture or leakage. Inspection techniques are needed in order to accurately estimate the minimum remaining wall thickness. A common approach is to perform a scan across the surface with an ultrasonic point thickness gauge, determining the wall thickness at each location from the time of flight of a bulk wave reflected off the back wall, assuming a known wave velocity. However, this is time consuming and suffers from the requirement to have access to the surface at all points; given that corrosion is often significant at supports, where water can collect and remain in contact with the surface, this access requirement presents a significant challenge.

Guided waves have long been developed for a number of applications including the inspection of pipes [1], rail [2], aircraft structures [3], adhesive joints [4] and even human bones [5]. They have the ability to inspect large areas from a single location, and technologies exploiting them for the long range inspection of pipes are now available as commercial products and sold worldwide [6]. However, while such techniques can successfully detect and locate damage within pipes, they are unable to accurately size defects, because the amplitude of the reflected guided wave, which can indicate severity, is dependent on several parameters beyond just defect depth, including the axial and circumferential extent. By contrast, guided wave tomography uses an array of transmitters and receivers at the boundary of the area of interest to collect data, and uses imaging approaches to reconstruct a thickness map from this. This approach was originally investigated by Jansen & Hutchins [7]. Other authors [8,9] have also undertaken extensive work to image damage using guided waves, although in the majority of cases the approaches do not quantify a physical parameter such as thickness.

Guided wave tomography is an inversion problem, to determine a map of wall thickness throughout the region of interest from a given set of guided ultrasonic wave measurements. The success of addressing this problem depends both on the type of guided waves (i.e. mode, frequency) used to interrogate the defect and the subsequent inversion technique. By far, the most common approach is to use Lamb modes, operating at a frequency at which they are dispersive. The dispersion curves, illustrated in figure 1, demonstrate how wave velocity is a function of the product of frequency and thickness, so by fixing frequency, velocity can be determined and then used to quantify thickness. Velocity estimation can be achieved by a number of techniques, most commonly ray tomography [10–13], with alternatives including diffraction tomography [14,15] or hybrid techniques [16]; these velocity inversion approaches were compared for guided wave tomography applications in [17], showing how the hybrid technique performed the best across a wide range of defects.

The widely used underlying assumption in guided wave tomography is that the thickness at every location maps to the corresponding velocity given by the dispersion curves. The dispersion curves are derived for an infinite plate of constant thickness, yet in guided wave tomography the interest lies specifically in cases where the thickness varies. Accordingly, Huthwaite [17] showed that this assumption was only valid where the thickness variation occurred over lengths of around 2 wavelengths or more for a 0.5 MHz mm excitation frequency with the A0 mode. This could be improved by an inversion approach which more accurately captures the guided wave scattering from thickness variations, rather than assuming it is the same as scattering from velocity variations. Full wave inversion, where a full numerical model of the guided waves is matched to the measured data, could offer a solution, but is slow and impractical for all but the smallest defects. Mindlin plate theory, investigated by Rose and colleagues [18–20], could be a solution, but has so far only been formulated for weakly scattering defects due to the use of the Born approximation. A completely different assumption to establish thickness is to use the cut-off properties of guided wave modes [21], but this has only been investigated for inspecting between two points rather than imaging an area.

Despite the strong potential of the velocity-to-thickness mapping approach, the understanding of the physical wave scattering from defects, and the consequences for guided wave tomography, is still limited. This paper aims to address this, firstly through a comprehensive study of reconstruction accuracy for the two fundamental modes, S0 and A0, across all appropriate frequencies, and then by measuring the arising error and investigating its source in terms of the guided wave scattering problem.

Following background theory in §2, the testing procedure is introduced in §3 and the results presented across both modes, for a wide range of different frequencies, using different defect widths and depths in §4. The role of effective density in these results is then discussed and assessed with some simple test cases in §5. Finally, in §6, an initial solution to correct for the density is introduced and its performance assessed for both simple and complex defects. The test data used in the paper is released to the wider research community as a benchmarking dataset to enable other researchers to directly compare their results, with details in §7.

## 2. Guided wave theory

Lamb waves [22] are guided waves travelling within an elastic plate. The plate had thickness *T*, set to 10 mm throughout the examples in this paper. To keep the results as general as possible, dimensions were normalized relative to *T*, so the results will be valid for any nominal plate thickness by varying frequency. The top and bottom plate boundaries were sound soft, i.e. behaved as if they were in a vacuum. In practice, air will surround the plate, which will behave in a very similar way to the vacuum as the impedance of air is so much lower than that of the steel plate. Figure 1*a*,*b* presents the dispersion curves for Lamb waves in a steel plate, calculated using the Disperse package [23,24], along with visualizations of the mode shapes for the two fundamental modes, A0 and S0 in figure 1*c*. These were derived using the assumption that the plate is of infinite extent and uniform thickness, and present the phase and group velocities, respectively, as the product of plate thickness *T* and frequency. The focus on this paper is on the fundamental A0 and S0 modes. Higher order modes can be used for guided wave tomography [25], but exciting and measuring pure modes, which is recognized as important for tomography, is challenging when higher order modes are present.

In this paper, mode conversion is neglected, in common with the majority of work on guided wave tomography thickness mapping. Pure A0 or S0 at frequencies below 2 MHz mm will mode convert into a combination of A0, S0 and SH0 when interacting with a defect, and the maximum amount of mode conversion will occur when encountering a very sharp 50% deep defect. In practice, the majority of defects are not like this, typically being a lot smoother, so mode conversion is likely to be minimal. The use of appropriate transduction to excite, and, by reciprocity, receive, pure modes also minimizes the effect of mode conversion; for interference from mode conversion to be a problem the waves must be mode converted twice (away from then back to the mode of interest) to be detectable. In practice, the most common effect would be a small energy loss from the mode of interest, manifesting itself as a localized attenuation.

## 3. Testing approach

To analyse the performance of guided wave tomography, imaging with the A0 and S0 modes across a range of frequencies and a number of defect widths and depths was performed, and the errors were quantified.

### (a) Defects

For the purposes of this initial study, simple axisymmetric Hann-shaped defects were used, the thickness maps of which were described by
*T* represents the background thickness, *D* gives the defect depth at its deepest point (the centre), *W* is the width and **r** represents the position on the plate. This is illustrated in figure 2*a*,*b*. While these defects are significantly simpler than true corrosion profiles, their simplicity greatly aids the analysis of the problem. It is also recognized that in previous studies (e.g. [17]) more complex defects have been found to behave as a superposition of these simpler defects. A more complex example will be investigated later in the paper. The defect depths were varied from 10% to 80% of wall thickness in 10% increments, and the defect width was set to 24, 16, 10, 8, 7, 6, 5, 4, 3, 2, 1.5 and 1 times the plate thickness; this totalled 96 defects.

### (b) Configuration

The configuration considered in this paper is shown in figure 2*a*. A circular array consisting of 240 transducers, and radius 40T was placed on the plate. These transducers can excite pure A0 or S0 across the full frequency ranges investigated (illustrated in figure 1*a*); clearly this is only possible for the finite-element (FE) simulations and, in practice, specific designs will be required to excite each mode at a particular frequency.

The defect was positioned at the centre of the transducer ring. Rather than acquiring 240×240 measurements for all the different send–receive combinations, requiring 240 simulations to model every source, only a single simulation needed to be performed for the case of a simple axisymmetric defect considered initially, and half of the measurements taken from the array recorded, as illustrated in figure 2*a*. The full dataset was then generated by replicating these signals.

### (c) Modelling

In order to capture the full complexity of guided wave scattering from corrosion defects, a full three-dimensional elastodynamic FE method was used; this does not make simplifications about the physics (beyond the FE discretization itself, which can be addressed through mesh refinement) and has previously been validated against experimental data for guided wave tomography problems in [17]. To maximize speed, the Pogo software package [26] (available from www.pogo-fea.com), which runs the explicit time-domain FE method on a set of graphics cards, was used.

The element size was chosen to be 1 mm (*T*/10) in each dimension, with eight-noded linear hexahedral elements used to model the plate, using the approach in [17]. As illustrated in figure 2*c*, to model the thickness variation associated with the defects, the elements were squeezed in the out-of-plane (*z*) dimension. In total, 15 million (1171×1171×11) nodes were used to model the plate, and the model was split between four Nvidia GTX Titan cards of 6 GB memory each. The time step was set to 3×10^{−8} *s*, which was smaller than ideal for a uniform plate because some elements had been reduced in dimension due to the plate squeezing, so the step had to be reduced in order to maintain stability of the FE method; 59 761 time increments were needed, and each simulation took around 50 min to run. To avoid reflections from the truncation of the mesh, some form of absorbing boundary was needed. For these purposes, a standard ALID (absorbing using increasing damping) boundary was used [27], although other, more optimized approaches are possible [28].

To excite the respective modes, forces were applied to a column of nodes through the thickness of the plate in the *z*-direction. To excite A0, these had uniform amplitude, but for S0, the force *F*(*z*)=2*z*/*T* was used, where *z*=0 corresponds to the central plane of the plate, to excite the symmetric mode; these are illustrated in figure 2*d*,*e*, respectively. It is noted that throughout the bandwidth of interest, only the fundamental modes are present, so pure antisymmetric or symmetric excitations (as described above) can only excite A0 and S0, respectively, provided any measurements are taken sufficiently far from the source to avoid non-propagating modes. In reception, to measure A0, the *z*-component at the central plane was measured. For S0, the difference between the *z*-components at the top and bottom surface was used.

A broadband time signal was applied so that all of the frequencies could be simulated with a single time domain simulation; this is achieved through a chirp excitation (e.g. [29]). The chirp, described by
*f*_{0} is the lower frequency, *B* is the bandwidth (i.e. difference between the highest frequency and lowest), *t*_{0} is the starting time of the signal, and *t*_{tot} is the desired signal duration. For the purposes of this study, *f*_{0}=30 kHz, *B*=170 kHz and *t*_{tot}=10/*f*_{0}, for both A0 and S0.

### (d) Signal processing

While the chirp excitation allows a very broadband signal to be excited, it requires some additional processing in order for ‘standard’ toneburst signals to be extracted. This can be achieved through a deconvolution approach, i.e.
*s*_{tb} corresponds to the desired toneburst input signal and *m*_{ch} is the measured signal from the FE data with the chirp excitation. The inversion algorithms required both arrival times and frequency domain values from the toneburst data.

As the guided waves are dispersive, distorting the wavepacket, determining a reliable measure of wave travel time is challenging. To help this, the guided waves are dispersion corrected to minimize their extent in time; for the A0 mode, this is achieved assuming propagation through a uniform plate, via the approach in [30]. For S0, this is less straightforward because the dispersion caused by the thickness variations of the defect is significant. Therefore, an additional correction was needed, and this approach is discussed in appendix B.

Frequency domain data at each frequency was obtained by performing a fast Fourier transform on the time traces. This can be undertaken on the chirp excitation to avoid the need to perform a separate Fourier transform. Following the approach of [16], these values were calibrated such that if there was no defect present, the values would be equal to the wavefield generated by a unit point source, as given by the Green's function
*r* is the distance from the source to the receiver. For the simulated data it is straightforward to do this by comparing the data to that of an incident field simulation, where there is no defect present.

### (e) Inversion

Following the standard guided wave tomography approach, a velocity map must be produced initially, prior to being converted to thickness. To do this, HARBUT (the Hybrid Algorithm for Robust Breast Ultrasound Tomography) [31] was used. This algorithm has been shown to produce accurate velocity maps across a wide range of different scatterer contrasts and sizes, down to a resolution limit of λ/2. It should be stressed that the goal of this paper is not to investigate the performance of this velocity inversion approach itself (this was investigated in [17]), but it is instead focused on the accuracy of the assumption that velocity maps to thickness. As shown in [17], this assumption results in a much poorer resolution limit (appox. 1.5 to 2λ) than that of the algorithm, so therefore the results presented here, using HARBUT, will be general for any velocity reconstruction approach whose resolution limit exceeds that caused by the thickness-to-velocity approximation.

The HARBUT algorithm combines the complementary strengths of bent ray tomography (BRT) and diffraction tomography (DT). BRT produces a reconstruction of velocity by fitting a ray model to the arrival times extracted from the measurements. An iterative approach is used, with the ray paths bending at each iteration to account for refraction. BRT is known to have a poor resolution limit [17,32,33]. DT is much higher resolution, but due to its reliance on the Born approximation [34], it is only suitable for weakly scattering objects, i.e. objects of low contrast and extent. HARBUT uses a BRT image as a background to reduce the residual contrast which needs to be imaged by DT. As discussed in [16], it is possible to improve the reconstruction by iteration, where the initial HARBUT image is used as a new background for additional accuracy improvements. This iterative version was used in this paper.

## 4. Parametric sweep results

Initially, each parameter was varied one at a time, while fixing others. Frequency variation was considered first of all, for a 40% deep defect, with width fixed to 8T, shown in figure 3. Figure 3*a* shows the A0 performance; as frequency increases, the accuracy improves so that the shape matches the true profile more accurately. The behaviour here is strongly indicative of an improvement in resolution as frequency increases; this would reflect the majority of imaging systems where resolution is proportional to wavelength, which reduces as frequency increases. The gain from increasing frequency seems to plateau at high frequencies, without the depth ever reaching the true value, suggesting that there is still a small residual which cannot be accounted for purely by considering resolution. Frequencies much higher than 1.2 MHz mm are challenging to invert because the group velocity becomes non-monotonic around this point. However, the presence of the plateau gives a strong indication that no further benefit in accuracy would be expected. For S0, the non-monotonic behaviour has necessitated the frequencies being split for clarity. In figure 3*b*, the lower frequencies are plotted. This time there is an overestimate in wall loss, and this overestimate reduces as frequency increases. That the error reduces as frequency increases is consistent with A0, but this time the primary cause cannot be resolution because the blurring caused by a poor resolution limit would cause the wall loss to be underestimated, not overestimated as seen here. Sensitivity is an issue with the S0 mode [35]; as its dispersion curve is very shallow at lower frequencies, a small error in velocity will result in a large error in reconstructed thickness. At higher frequencies, the gradient of the dispersion curve is greater so this is less of a problem. This is clearly visible in figure 3*b*, with the error reducing with increasing frequency. Figure 3*c* appears to contradict this, however. As frequency increases further beyond 1.8 MHz mm, the error becomes larger, despite figure 1*a* showing that the gradient increases more at higher frequencies. This is a reflection that the error is not just a constant offset, but is likely to be proportional to velocity contrast in some form. At higher frequencies, therefore, the higher velocity contrast caused by the steeper curve causes a higher error, and because the error will be projected along the flatter section of the curve, it will have a disproportionate effect on the resulting thickness. This analysis would suggest that an operating point of approximately 1.7–1.8 MHz mm would give the best results for S0 for this particular defect size and shape.

Figure 4 presents the results for a 40% deep defect as the width was varied, fixing frequency at 0.5 MHz mm for A0 and 1.8 MHz mm for S0. Considering A0 initially (figure 4*a*), there is a slight underestimation error visible with the largest defect, matching the underestimate of figure 3*a*. As the width reduces, the error increases; this is likely to be related to the limited resolution as discussed earlier. For the S0 plots (figure 4*b*) once again there is an initial offset for the widest case (16T), but this shows an overestimate of wall loss, which matches what was seen in figure 3*b*,*c*. As the width reduces, this overestimate error increases, to a peak of around 12% overshoot at a width of 8T. From this point the reconstructed depth reduces monotonically with reducing width, as for the A0 case, reflecting the loss of resolution. It is clear, however, that the increase in error from 16T to 8T cannot be caused by loss in resolution, because limited resolution will cause blurring and hence underestimation, not the overestimation seen. There is clearly another contribution to the error, and this will be discussed later in the paper.

The variation of imaging performance across depth is shown in figure 5, when width is 8T. For A0 at 0.5 MHz mm, the results are shown in figure 5*a*, illustrating how the reconstruction depth increases as the true depth increases. The linearity of this behaviour is clarified in figure 5*c*, where the reconstructions are represented as wall loss and expressed as a fraction of the true wall loss. From this, it is clear that the error seen is proportional to the depth; the lines are difficult to distinguish, and the error appears to be approximately 30% of the defect depth. For S0, shown in figure 5*b*,*d*, the results appear less linear, with the error fraction increasing as wall loss increases, until the remnant thickness is clipped at 0%. This nonlinear behaviour is a result of the dispersion curve being highly nonlinear around frequency-thickness values corresponding to these defects.

The previous analysis has enabled a single parameter to be varied while fixing all others. To generalize this, an error measure was taken. In [17], analysis was performed with both the error in the maximum depth and a weighted RMS error across the image. It was found that these two metrics showed the same trends, so for simplicity, just the former is used here. Error is expressed as a fraction of the nominal wall thickness, *T*, and the ‘signed’ values are recorded such that it is clear if it is an overestimate (positive) or underestimate (negative) in wall loss. Figure 6 shows the results from this, plotting the error values for A0 and S0, across a range of frequencies and depths, for four different defect widths.

In figure 6*a* (16T wide, A0), there is good performance for the majority of depths across most of the frequencies. However, at very large depths (70% and 80%) there are large errors; these occur because of the very high contrast of the defect combined with its width, making it fall very far outside the limits of the Born approximation underlying DT. In these cases, BRT does not compensate enough to enable the ‘global minimum’ to be converged to. It is noted, however, that the lack of convergence can clearly be seen by large jumps between iterations, so there is no danger of incorrectly identifying the resulting image as being correct; in all the results presented in figure 6 this is the only error which has been caused by inaccuracies in the velocity reconstruction algorithm. A common approach to address these convergence issues is to iterate from low frequencies to high frequencies, which has not been exploited here; the results are for single frequencies. Moving up the *y*-axis to the lowest frequencies, the error is clearly increasing, matching the trend observed in figure 3*a*. This monotonic trend is visible across all the defect widths with A0 (figure 6*a*,*c*,*e*,*g*), and is primarily caused by the lack of resolution blurring the reconstruction. The other trends seen for A0 in figures 3, 4 and 5 can be identified from this figure; the error increases both as width reduces and depth increases.

In the previous analysis for S0, an optimal frequency for the specific defect (40% wall loss, 8T wide) was identified to be around 1.7–1.8 MHz mm. From figure 6*b*,*d*,*f*, S0 reconstructions at widths of 16T, 8T and 4T, respectively, this value holds as producing the best results across all depths. The tendency for the reconstruction depth to saturate the wall thickness with S0 at frequencies above and below around 1.7–1.8MHz mm is clear from figure 6*b*,*d* particularly, with the constant (as a function of frequency) values corresponding to a reconstruction of 100% wall loss, and hence an error of *T*−*D*. In figure 6*f*,*h*, underestimates are visible, as the defect widths are smaller and the limited resolution effects appear in the reconstruction.

One of the clearest outcomes is that A0 underestimates wall loss, while S0 overestimates, certainly for the wider defects where resolution effects are not an issue. This systematic difference is potentially useful in that it could be exploited, with A0 and S0 used together to define a range over which the true defect depth may be. However, by analysing and understanding this it is possible to correct for this term, therefore, improving the accuracy of reconstruction from a single mode. As the next section will show, this term is related to scattering from a ‘density type’ term, where the density is proportional to the thickness of the plate.

## 5. Density in guided wave scattering

When trying to analyse the errors observed previously, it should be noted that the underlying assumption, that thickness maps to velocity according to the dispersion relationship, must hold if there are no thickness changes. Therefore, for a sufficiently large region of constant thickness, no error should be visible, assuming that the dispersion relationships are valid, the test data have been produced accurately and the velocity inversion has been successful. This leads to the conclusion that any errors, such as the overestimate for S0 and underestimate for A0, must be associated with localized variations in the thickness.

The object function is the underlying quantity reconstructed by an imaging algorithm. In the case of diffraction-tomography-type algorithms, the object function can be defined, including density, as [36]
*k*_{u} and *c*_{u} correspond to the constant (uniform) background values for wavenumber and velocity, respectively, *c* is the local velocity and *ρ* is the density. The use of velocity within this acoustic model for guided wave tomography is very well established, being mapped to thickness via the dispersion curves. A logical definition of the effective density for guided waves will be to assume that it is proportional to the mass per unit area of plate, i.e. proportional to thickness, assuming that there are no variations in material density. It is noted that the form of the equation means that any constant of proportionality in the density term will automatically be cancelled, so this does not need to be considered.

Importantly, this density term would fulfil the criteria outlined above for the error seen; for regions of constant thickness, the density term of equation (5.1) will be zero due to the Laplacian (∇^{2}) operator. Also the density term would be similar for A0 and S0, certainly with the same sign, whereas the velocity contrast will have opposite signs, so therefore, the density term is likely to lead to an overestimate with one mode and an underestimate in the other.

The simplest case for testing this is an axisymmetric flat-bottomed-hole with a taper down from the background area to the flat region; this is illustrated in figure 7*a*,*b*. In this model, the flat bottom itself is a constant thickness region, so should have no error associated with it. As described previously in §3, an FE model of this defect is run, using A0 at 0.5 MHz mm, to generate a set of test data, and a velocity map is produced from this dataset, using HARBUT, which is then converted to a thickness map. Following the convention of [17], the test dataset produced from the full three-dimensional elastic guided wave model will be referred to as ‘Type 1’, as will the corresponding images produced from this dataset.

Two new datasets were then produced for comparison from a finite difference (FD) model. This is a two-dimensional acoustic FD code, running in the frequency domain. This model was an 800×800 grid of nodes, with the nodes separated by 1 mm. This was formulated to discretize the acoustic wave equation including both velocity and density variations. Full details of this model are given in [17]. ‘Type 2’ data was generated from the FD model where density was constant and velocity varied to match the dispersion curve value at the local thickness. A third dataset was introduced, ‘Type 3’, identical to Type 2, except with the inclusion of density variation, proportional to thickness. Reconstructions of thickness are produced from Type 2 and Type 3 data in the same ways as Type 1.

To analyse the results, differences between the reconstructions were taken. As the generation of Type 2 data matched the inversion approach (a pure velocity model), this was taken as the ‘baseline’ measurement, and the reconstructions of figure 7*c*,*d* are given as reconstructions of Type 1–Type 2 and Type 3–Type 2, respectively. Note that in all cases, to remove the effect of resolution, all images were filtered to remove spatial frequency components above 1.6λ, reflecting the resolution limit discussed in [17], and similarly any small DC offset is removed (this occurs due to the singularity in the beamforming to diffraction tomography filter [37]), prior to the image being windowed to remove any artefacts which appear around the transducers. Figure 7*e* plots cross sections through the centre of the defect for these two plots, along with the true wall loss.

The Type 3–Type 2 image, plotted in figure 7*d*, shows what effect the density term of equation (5.1), −*ρ*^{1/2}∇^{2}*ρ*^{−1/2}, would have on a thickness map. In effect, this could have been generated directly by calculating the object function caused by the density term in isolation, calculating what velocity map would correspond to this object function, then mapping this to thickness. Owing to the Laplacian operator, there are spikes seen when the gradient of thickness/density is discontinuous; this is most clearly seen by the comparison with the wall loss plot in figure 7*e*. While the Laplacian for such a jump mathematically would be a sharp discontinuity, the images presented are limited by resolution which blurs the spikes.

The Type 1–Type 2 reconstruction, plotted in figure 7*c*, shows a very similar behaviour to that of the density residual of Type 3–Type 2. Clearly, there are two concentric rings in the same locations shown in the images, and the corresponding spikes in the cross sections are in the same positions as shown in figure 7*e*. However, both the sign and the amplitude are different, which will be analysed below.

### (a) Sign

Figure 7 showed that the sign of the residual density-type term reconstructed in the image differed between the guided wave scattering model and the acoustic model. To conduct further analysis, very low frequency S0 is considered, where the dispersion curve is approximately constant. This has the benefit that the guided wave scattering is almost entirely caused by the effective density change, rather than velocity variation. If the sign change is visible with a guided wave at this frequency, it would help provide an explanation for why a sign change appears elsewhere both with this mode and others.

The set-up considered is very similar to the case considered before, except with an excitation frequency of 0.3 MHz mm for the S0 wave, which resulted in a much larger wavelength. Owing to this, the model was scaled such that the array radius was 1.8 m, and the defect outer and inner diameters were 1.45 m and 0.55 m, respectively. The elements were resized to 4×4 mm^{2} in the plane, and 1.25 mm tall; while this resulted in elongated elements, it was important to keep a number of elements through the thickness to capture the full mode shape. For comparison, a Type 3 model (i.e. acoustic including density variations) was run, along with a two-dimensional FE plane stress elastic model, where the material density of each element is proportional to thickness, with the Young's modulus also being reduced by the same proportion to avoid velocity changes.

Figure 8*a* is the thickness map used, also including the circular array; transducer 1 was the source. A complete set of traces for this source are presented in figure 8*b*,*c*,*d*. Four separate wavepackets are visible away from the through-transmission point, marked a–d. These correspond to the wave interacting with four boundaries, as illustrated in figure 8*a* between the source and receiver 75. The fully three-dimensional guided wave model data is in figure 8*b*, and the important feature to note is that the wavepackets are not continuous from reflection through to transmission; each one passes through a ‘zero’ (marked). By contrast, in the results from the two-dimensional acoustic model of figure 8*c*, the waves are continuous. This is a reflection on the nature of acoustic scattering; typically, a point scatterer will behave as a monopole, while the indication of figure 8*b* is that the guided wave scattering is closer to a dipole. This behaviour has been previously recognized by Rose *et al.* [19] for guided wave scattering. Importantly, it is also visible in two-dimensional elastic scattering, shown in figure 8*d*, where once again there are clear zeros visible in the wavepackets at particular angles.

Figure 8*e*,*f* compare the time traces for the scattered field at transducers 1 and 75, respectively. For transducer 1, in direct reflection, it is clear that all three time traces match each other well, with no phase shifts between signals, indicating that in reflection, a thickness change behaves in the same way as the corresponding density change. For transducer 75, as shown in figure 8*f*, while the guided wave and elastic signals match each other, they deviate significantly from the acoustic signal. This deviation depends on the specific path that has been taken. As shown in figure 8*a*, the wave paths a–d each correspond to different deviation angles, from a which is closest to transmission with very little deviation, to d which is closer to reflection. For the ‘reflection’ wavepacket d, it is clear that while the amplitude of the guided wave is reduced compared with the acoustic wave, there is no phase shift; it is similar for wavepacket c, although the amplitude here is further reduced. For wavepackets a and b, which can be considered to fall into the ‘transmission’ subset, along with an amplitude drop there is also a sign change. This sign change is an indication of the dipolar scattering behaviour discussed earlier for guided waves compared with the monopole scattering for acoustic data. The transition in sign will be associated with a zero as seen in figure 8*b*.

In guided wave tomography, as discussed in [17], reconstructions are performed with a subset of transmission components close to through-transmission. The component scattered by density variations within this transmission subset will have a sign change when compared with the equivalent acoustic model, as shown from the analysis above. However, the acoustic model has been used to derive the density term *ρ*^{1/2}∇^{2}*ρ*^{−1/2}, and therefore, a sign change is needed for it to be used when accounting for density scattering with guided wave data. It is this that explains the sign difference between figure 7*c*,*d*.

### (b) Amplitude determination

Having established that the dipolar scattering behaviour gives a reason why the sign of the density term is inverted in the guided wave reconstructions relative to the acoustic equivalent, there is still the issue that the absolute amplitude of the reconstructed density component in figure 7 differs from the value suggested by the acoustic model. As shown in figure 8, the scattering behaviour for the guided wave case was shown to vary in sign and amplitude depending on scattering angle, compared with the acoustic case. Here, the amplitude of the guided wave field scattered by density will be taken as a constant (with scattering angle) factor, *Q*, multiplied by the equivalent acoustic wave field. It is assumed that the true variation of this factor through the relatively narrow subset of transmission data used can be neglected. By comparison, it should be noted that *Q*=1 would be expected for the reflection component, as these reflection components were shown to match in the previous section.

In this section, the aim is to determine this factor, relating the transmission scattering behaviour of guided waves from density variation to that of acoustic waves. The options to determining this ‘effective density’ factor include (1) analytically analysing the problem, e.g. through Mindlin theory or similar approaches, to determining the scattered response, (2) comparing the scattered field in a similar manner to figure 8 or (3) imaging the scattered field and extracting the density component from that, similar to what is shown in figure 7. To build up a general picture, (1) is likely to be challenging; while models exist which can be used to do this, they are generally limited in applicability and would struggle to provide a meaningful data for both A0 and S0 across the full frequency spectrum. Extracting the relatively small density component in isolation from a measured signal will be challenging, so (2) is similarly eliminated.

The remaining option is therefore (3), comparing the density components of the images produced. The density term is relatively small (yet as shown in figure 6, important in the accuracy of the reconstructions), and hence measuring it is dependent on its easy separation from the velocity component. The axisymmetric Hann shape used earlier suffers in this regard; the largest value of ∇^{2}*ρ*^{−1/2} will occur at the central point where the velocity contrast is also the largest, making the separation challenging. Instead, a large flat bottomed hole with a taper is once again proposed, of the type investigated in figure 7. In this, the shape of ∇^{2}*ρ*^{−1/2} will be very distinct from the shape of *c* in the reconstructed image, so this should aid in their separation.

#### (i) Model

A full elastic FE model was run in the same manner as before, for both modes and with a chirp excitation to enable multiple frequencies to be extracted. It was necessary to make the flat bottomed hole as wide as possible because the poor resolution limit for long wavelengths (particularly, S0 at low frequencies) can otherwise hinder the separation of the velocity and density components; for this reason, the mean width was taken as 600 mm with a taper of 250 mm. The depth was set to 1 mm; while a deep defect may help to maximize the density component above any artefacts present in the image, such defects can lead to greater errors in the resulting image because of nonlinearities in the inversion approach. One of the largest sources of errors with the HARBUT algorithm is that errors occur around the transducer locations; while these can be gated, it is necessary to make the array radius as wide as possible. The modelled radius was 800mm, however, the field was then analytically propagated outwards to give a wider radius of 2400 mm; the details on how this was done are given in appendix A.

The plate was 10 mm thick as before, and in order to accommodate the large array and the 195 mm absorbing boundaries needed to absorb the broad bandwidth of waves, the width was 2 m. Meshing with 1 mm elements created 44 million nodes, requiring about 70 GB of GPU memory. For this reason, the simulation was performed using Pogo [26] on a system with eight dual K80 Nvidia GPUs. The 48 934 time steps took around 1 h to complete using six GPUs for each simulation.

#### (ii) Component separation

Following the definition of equation (5.1), the object function should be separable into two components
*O*=*O*_{vel}+*QO*_{den}, where *Q* is expected to depend on mode, frequency and local thickness, and is to be determined.

An automatic fitting approach can be used to achieve this, but this presents challenges because the presence of artefacts in the reconstruction can mask the true amplitude. In recognition of this, therefore, the coefficients are fitted by manually adjusting the multiplier on *O*_{den} until it matches the residual *O*−*O*_{vel}. The results are given for A0 and S0 in tables 1 and 2, respectively.

It should be recognized that the accuracy of the values is limited by the presence of artefacts which occur when considering these small components. Within the expected accuracy, the results suggest a lack of variation with frequency. In the case of S0, the artefacts became too severe above 1.7 MHz mm to enable the components to be separated, however, given the results in §4 suggesting the optimal operating point is around this limit, this should not pose a significant problem. For the ranges given, the multipliers for A0 and S0 for shallow defects are taken as *Q*=−1.9 and *Q*=−0.33, respectively. As discussed, in reflection *Q*=1 is expected, so these different values for *Q* in transmission suggest scattering patterns with different shapes for the two modes.

#### (iii) Variation with defect depth

The relative invariance of the scaling parameter with frequency is a convenient simplification which can be exploited when considering varying depth. As a thinner section of plate can be considered to behave in the same way as a thicker region at a lower frequency, initially it might indicate that the scaling factor will remain constant as depth increases. However, the Laplacian operator will need to be adjusted; the amount of scattering will be determined by the rate of thickness change per wavelength, so the derivatives forming the ∇^{2} term in the equation must be adjusted accordingly. Based on this, the object function equation can be expressed as
*c*/*f*

## 6. Density correction

Having established a model describing how density will affect imaging within guided wave tomography, this can now be exploited to improve the imaging performance. Given that neglecting density still enables a reasonably good image to be produced in the majority of cases, as shown earlier, the object function can be assumed to purely be a function of velocity, with the initial estimate being referred to as *c*_{1}, i.e.
*t*_{1}
*f* is the frequency, *c*=*D*( *ft*) is the dispersion curve giving a velocity for a particular frequency–thickness product, and *D*^{−1} is its inverse. As the absolute value of *ρ* is immaterial because any constant multiplier will cancel out in equation (5.5), here it is simply normalized such that it is 1 when thickness is equal to *T* to give the effective density term

### (a) Application to simple defects

The approach outlined in the previous section has been applied to the simple Hann defects studied before, and some results are presented for a number of widths and depths, in figure 9, with the A0 mode at 0.5 MHz mm.

The results show a significant reduction in error across the full range of defects studied. It is noted that as the defect width reduces to 8T, the accuracy of the corrected case notably reduces; this is caused by the limited resolution discussed earlier. In figure 9*d*, showing how the error varies with depth, the correction improves the accuracy significantly. There was a reduction in accuracy for the deepest cases, but this is to be expected due to the high contrast providing a more significant challenge to the inversion algorithms.

To confirm that this approach works with experimental data, the data from [17] were used to produce a pair of thickness maps, both with standard HARBUT and with the density correction. The resulting cross sections are given in figure 10. These confirm the results from the FE models before, showing that when the algorithm is not resolution limited, density can be corrected to significantly reduce the maximum error.

### (b) Application to complex defects

The same principle can be applied to more complex defects. A laser scan of a real corrosion patch was used to generate a FE model, as described in [17]. The data were used to generate both a standard HARBUT image, and a HARBUT image with density correction. Figure 11 shows these results. Figure 11*a* presents the true thickness map, which was used in the FE model to generate data. For reference, as the inversion algorithm is resolution limited, a filtered version of the true map is given in figure 11*b*; this is filtered to remove any components above the 2λ resolution limit. The standard HARBUT reconstruction is shown in figure 11*c*, with the density-corrected version in figure 11*d*. The cross sections are plotted in figure 11*e*.

Applying density correction to HARBUT has significantly improved the reconstruction accuracy, with the values getting closer to the true filtered version, the theoretically best achievable result. There do seem to be some spurious values, however, such as the peak which appears in the middle of the cross section in figure 11*e*. It may be that for more complicated defects, the inversion of density and velocity terms in this manner represents a poorly conditioned problem. Other possible approaches to help the separation could be to exploit multiple frequencies; as the two terms differ in their behaviour as frequency varies this may provide an approach to separate them (e.g. [38]).

## 7. Benchmarking dataset

It is recognized that there are many inversion approaches which have been proposed by researchers to produce thickness maps from guided wave tomography data. This paper has used HARBUT as a vehicle for investigating the density term, which has been compared previously with two other commonly used algorithms [17]. Comparing techniques between different researchers is challenging because of the use of different defects and methods to generate the data. In an attempt to address this, the dataset used to generate the results in this paper will be made available for benchmarking purposes. The repository of data can be found at http://dx.doi.org/10.5281/zenodo.44626. As the entire dataset in its raw form is around 17 GB, it has been necessary to select specific cases; for this the defect width/depth pairs (both in millimetres for the 10 mm plate) are 3/100, 4/30, 4/80 and 5/240. This provides the data for the cross sections in figure 9 in addition to a narrow case (30 mm wide) in order to test resolution. In each case, data for both A0 and S0 modes are provided. Full data for the chirp input is provided, along with a Matlab function to enable the extraction of wavepackets at specific frequencies. Additions or suggestions from other researchers would be welcome.

## 8. Conclusion

The simulation and imaging of a range of defect widths and depths across both the A0 and S0 modes across a range of frequencies has shown a number of trends. For A0, an underestimate of wall loss is visible in the reconstruction. Resolution, and hence accuracy, improves with increasing frequency, and the error was shown to remain proportional to defect depth. In contrast to the behaviour of A0, S0 typically overestimates wall loss. The dispersion curve characteristic is such that there is a very clear point around 1.7–1.8 MHz mm where the accuracy is best, and away from that, large errors occur, suggesting that only a narrow range of frequencies is suitable.

A key outcome of this paper is to highlight that the reason for the A0 underestimating and S0 overestimating is the presence of an additional scattering term, related to the effective density. This effective density is proportional to the thickness, and scattering occurs from locations where the Laplacian (∇^{2}) of density is non-zero. While the standard acoustic model of density scattering needs to be adjusted for guided waves via the use of a scaling constant, it is possible to correct for this term in the reconstructions. As the error in many cases is almost entirely caused by neglecting the density, accounting for it has been shown to result in a near-elimination of error in many reconstructions, including when using experimental data and imaging complex defects.

It should be observed that the Laplacian will be largest when thickness changes rapidly, and this is likely to occur at the deepest point. This point is the most important to accurately reconstruct because it represents the minimum wall thickness and hence directly informs component life decisions. Therefore, correcting for the density term at this location is extremely important. The conclusions of this paper, highlighting the presence of the density term and demonstrating how it can be accounted for in a reconstruction, provide a step change in the accuracy of component life estimation as determined by guided wave tomography.

## Data accessibility

Data supporting this article can be downloaded from http://dx.doi.org/10.5281/zenodo.44626 as described in §7.

## Authors' contributions

The sole author undertook 100% of the work in this paper.

## Competing interests

The author has no competing interests.

## Funding

Funding for this work was provided via EPSRC grant no. EP/M020207/1.

## Acknowledgements

The author would like to thank Prof Peter Cawley and Prof Michael Lowe for their helpful discussions and critical reading of the manuscript.

## Appendix A. Analytical propagation

A scattered field (travelling outwards) can be expanded as [39]
*n*th Hankel function of the first kind and *a*_{n} represents a set of constants. The constants *a*_{n} typically drop to very low values quickly away from 0, so the series can be truncated. By sampling *ψ* at a constant radius at equidistant points around the full range 0≤*θ*<2*π*, it is clear that (8.1) becomes a Fourier transform. Therefore, it is straightforward to determine *a*_{n} from such a set of data. This can then be propagated to a different radius by changing *r*, and inversing the Fourier transform to give *ψ*. This can change the radius of the measurement locations. It is, however, possible to perform the same position alteration with the source locations, due to the principle of reciprocity.

## Appendix B. S0 arrival time extraction

Figure 12*a* shows an S0 wave signal, excited at a centre frequency of 1.75 MHz mm and having propagated 80T through a Hann defect 24T wide and 0.1T deep. Also plotted is an ‘incident’ signal which was propagated 80T without the defect being present, and the input signal. These are analytically generated by Fourier synthesis, simply summing the phase shift along the path at each frequency then performing an inverse Fourier transform. Clearly, the significant dispersion present acts as a hindrance to extracting the arrival time; the two signals are extremely elongated and the shape is likely to have a much greater influence than the true arrival time.

Dispersion can be minimized for in a number of ways (e.g. [30]); here an initial approach is to divide the signal at each frequency by *l* is the propagation distance and *k*_{0}( *f*) is the background wavenumber at the specific frequency. This corrects for propagation in the background medium, but does not account for thickness variations. The results are shown in figure 12*b*. While the dispersion is significantly reduced, there is still notable distortion which would mask the arrival time. This arises from the unknown variations in thickness associated with the defect which are not correctly accounted for with the uniform thickness assumption.

To achieve better results, it is therefore necessary to correct for the distortion caused by the defect itself. Neglecting reflection and scattering (assumptions which are consistent with bent ray tomography), the measured field at a particular frequency *f* can be expressed as
*k* along the ray path. The wavenumber is a function of frequency, *f* and thickness Th which depends on position. The background component (wavenumber *k*_{0}) can be separated from this as before
*c*, or its reciprocal, slowness, *s*
*T*, has been used. One further approximation is then used, that the slowness dispersion curve can be linearized, so that for any frequency there is a constant *K* where
*K* at each frequency. Based on this
*K* should be constant with frequency; to, in effect, force this to be the case, the frequency values can be mapped to new frequency values
*f*_{0} is the centre frequency, so
*ϕ*′ can then be inverse Fourier transformed to give an undistorted signal, within the linear approximation chosen. Figure 12*c* shows this result. Here the arrival time shift is very clear to identify from the signal. Because equation (8.8) normalizes against the centre frequency, the time shift will correspond to the phase velocity at this frequency.

- Received December 7, 2015.
- Accepted January 14, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.