## Abstract

Producing accurate thickness maps of corrosion damage is of great importance for assessing life in the petrochemical industry. Guided wave tomography provides a solution for this, by sending guided waves through the region of interest, then using tomographic imaging techniques to reconstruct the thickness map, importantly eliminating the need to take measurements at all points across the surface. However, to achieve accurate maps, the imaging algorithm must account for the way in which the guided waves interact with corrosion defects, and the complex scattering which occurs. Traditional approaches have exploited the dispersive nature of guided waves: a velocity map is produced from a dataset, then converted to thickness using the dispersion relationship. However, these relationships are derived for plates of constant thickness, which is not the case in the majority of defects, causing significant inaccuracies to exist in the images. This paper develops a more sophisticated inversion solution which accounts for the full-guided wave scattering, enabling more accurate images with resolution better than a wavelength, compared with two wavelengths previously. This is demonstrated with simulated and experimental data. The speed and stability of the algorithm in the presence of random noise and systematic errors is also demonstrated.

## 1. Introduction

Assessing the remaining life of assets is critical in the petrochemical industry, particularly in the presence of corrosion which reduces wall thickness of pipes and pressure vessels, increasing the chance of unexpected failure. Guided wave tomography [1] is a quantitative technique proposed to produce wall-thickness maps of pipes and plates; this enables the remaining thickness to be extracted and hence the life of the component to be determined, without requiring access to all points on the surface, making it a potentially valuable non-destructive testing technique. Lamb waves are excited in a plate or plate-like-structure (including pipe walls provided the wall is thin compared with the radius of curvature, meaning that curvature can be neglected [2]) from an array of transducers; these then interact with a region potentially containing wall thinning caused by corrosion, and the resulting waves are measured by another array. These measurements are subsequently processed to produce a thickness map.

Clearly, it is critical that these maps are an accurate representation of the true map to avoid the risks associated with overestimating remaining life or waste associated with underestimating it. This accuracy is highly dependent on the underlying assumptions used in the reconstruction. The primary assumption used reformulates the guided wave scattering problem as one of the acoustic waves being scattered by velocity [3]. As shown in the dispersion curves of figure 1, the guided wave velocity is a function of the frequency-thickness product, so by fixing frequency, wave velocity can be shown to be a function of thickness. This is a widely used assumption which has been used for 20 years in guided wave tomography [3–12] but relies on a key approximation. The dispersion curves are derived for infinite plates of constant thickness [13]; therefore, the assumption is made that any local thickness variations are sufficiently gradual that they can be neglected. As highlighted in [14], this manifests itself as a resolution limit in the reconstruction of around 2*λ*, effectively applying a low-pass filter to remove any higher spatial frequency components. This limit corresponds to 8*T*, where *T* is wall thickness at 0.5 MHz mm *A*_{0}, or 4*T* at the constant group velocity around 1.4 MHz mm *A*_{0}. Importantly, this limit is sufficient to cause significant inaccuracies with typical corrosion defects.

This lack of resolution when compared with the traditional limit of *λ*/2 can also be considered an effect of the different scattering behaviour observed between guided waves with varying thickness and acoustic waves with varying velocity. Figure 2 presents a schematic of this. Figure 2*a* defines the scattering angle, while figure 2*b* plots the scattering patterns for an acoustic and guided wave model. The acoustic case has an omnidirectional behaviour, whereas the guided wave model is more complex. From the diagram, it can be seen that the transmission components, which provide information about the low spatial frequency components of the image (e.g. [15]), match between the guided wave and acoustic models, enabling a low-resolution image to be produced. However, in other scattering directions (closer to reflection, corresponding to high spatial frequency components), there is a significant difference in scattering behaviour, which prevents higher resolution components from being extracted for imaging.

Significant focus has been given to the development of velocity inversion schemes, using this assumption that acoustic scattering matches guided wave scattering behaviour. Dominating this is the use of ray theory in the inversion. Typically, arrival times are extracted from the measured waves, then a ray model is fitted to these such that the modelled and measured arrival times match. This can use straight ray assumptions [1,5,16], although more recent approaches account for refraction by iteratively bending the ray paths [9,11]. Resolution improvements are available through the use of diffraction tomography for velocity inversion [8], although this is only applicable to weakly scattering defects, i.e. ones which are low contrast and small. The HARBUT approach [17] overcomes this by combining ray tomography with diffraction tomography to enable the benefits of diffraction tomography to be realized across a wider range of defects, and was applied to guided wave tomography in [10]. Full wave inversion, where an acoustic finite difference model is iteratively updated to match the measurements, has also been used for this problem [12], although the approach is still restricted to acoustic, rather than guided wave, scattering and, therefore, is limited in the same way as other velocity inversion approaches.

To improve accuracy beyond this, approaches that can more accurately capture the guided wave scattering are needed. A simple enhancement is to consider density in addition to velocity when inverting [18], which can improve accuracy, but this still fails to capture guided wave scattering fully. In a more sophisticated approach, Rose *et al.* developed Mindlin plate theory [19] for guided wave tomography [20,21]. Under this model, the behaviour of the plate is modelled using three field parameters: two rotations and the out-of-plane displacement. The scattering from guided waves can be captured by this model, and the paper illustrates scattering patterns and discusses how this behaviour can be used in an inversion under the Born approximation to achieve a quantitative map of thickness. However, the approach assumes small perturbations in local properties, i.e. small wall losses, and the use of the Born approximation severely limits the range of defects which can be modelled. Additionally, while Mindlin theory is undoubtedly an improvement on the guided-wave-to-velocity mapping assumption, it is a simplified model and, therefore, cannot capture the full scattering behaviour. Iterative solutions, such as the Monte Carlo Markov Chain, have been applied to this problem [22], but this was limited to very small scatterers due to the need to run the forward model many times. Other guided wave imaging approaches have been demonstrated, including the use of the linear sampling method [23], and techniques producing non-quantitative images [24].

This paper develops a new general model, fully capturing guided wave scattering, and develops an inversion technique based on this. The approach is based on the HARBUT algorithm [10,17], makes no restriction on the width or depth of the scatterers and generally provides results that have superior resolution and accuracy to the existing techniques. Section 2 discusses the theory, providing the background of the standard acoustic (i.e. non-guided-wave) scattering model which is then extended to enable scattering from guided waves to be captured. Section 3 then outlines approaches to invert the developed scattering model, with §4 presenting results from the algorithm. Section 5 discusses aspects such as stability and effects of density.

## 2. Scattering theory

### (a) Standard acoustic model

The analysis begins with the standard acoustic model, widely used for many quantitative imaging problems, from which a more general approach will be derived which better captures complex scattering. The acoustic wave equation is written in the frequency domain
*k*(**x**) is the local wavenumber at position **x** and *ϕ* represents the local acoustic scalar potential. This can be rearranged
*k*_{0} is the background wavenumber, and the object function *O*(**x**) has been introduced; this is normally expressed in terms of the velocity as
*c*_{0} represents the background velocity and *c*(**x**) the local velocity.

The acoustic free space Green’s function can be derived as the solution to
*ϕ*_{s} is derived in terms of the total field *ϕ* and the incident field *ϕ*_{0}.

Taking a point scatterer of scattering potential *ρ* at position **x**, the integral reduces to
*G*. It is well known that guided waves do not scatter like this [26,27]. As shown in [14], it is the discrepancy between this acoustic scattering and the full-guided wave scattering which limits the accuracy of the acoustic model and hence its resolution. To improve this, a more accurate scattering model is required.

### (b) Generalized scattering model

As discussed, the monopole scattering behaviour is unrealistic for capturing full-guided wave scattering. Therefore, higher orders are needed; these can be expressed as a sum of cosine terms in a variation on equation (2.7)
*θ* is the deviation angle between the illumination and reception direction and *M* is the highest order cosine term; for the omnidirectional acoustic case *M*=0 is sufficient. The values *C*_{m} represent the coefficients associated with each cosine term and hence describe the scattering behaviour. It is noted that the added directionality term could be included in a Lamb wave Green’s function; however, the terms will be kept separately throughout this paper and *G* will always be used to refer to the two-dimensional (2D) free space acoustic Green’s function.

Defects across a wide range of depths (including beyond 50% wall loss) are of interest so it is not guaranteed that the scattering pattern will remain constant with depth. Therefore, the model must be adjusted so that different scattering patterns can be expressed at different depths, which is achieved via a polynomial expansion of the object function in terms of the thickness parameter, *h*(**x**)=[*T*_{0}−*T*(**x**)]/*T*_{0}, with *T*_{0} and *T* being background and local thicknesses, respectively,
*C*_{mn} has been expanded to account for the *N* polynomial terms. This form of equation allows the scattering pattern to vary with depth. It is noted that, as scattering is not expected to originate from the background region where *h*=0, there is no constant term, so *n*=0 is not included. While this equation represents an approximation due to the truncation of both the polynomial (*n*) and cosine (*m*) series, only low orders of both are typically required to capture the scattering behaviour.

To predict the scattering behaviour under this model, the coefficients *C*_{mn} must be estimated. Rearranging the equation gives
*ϕ*_{s} and *P*_{mn}, a suitable forward model is needed. Traditionally, significant simplifying approximations were necessary to produce such data. However, the increasing computer power available makes running large sets of full finite-element (FE) simulations practical. The Pogo FE software package [28] solves elastodynamic problems via an explicit time domain scheme, and is used here to simulate the wave interaction with defects. This approach captures the full-guided wave scattering behaviour, enabling the field *ϕ* to be determined. By running a similar simulation without the scatterer present, the incident term *ϕ*_{0} can be determined and subtracted from *ϕ* to give the scattered field *ϕ*_{s} away from the scatterer, and hence the left-hand side of (2.12). The integral (2.13) is calculated using a sum across the domain **x** with a sufficiently small spacing (typically *λ*/15 is sufficient, depending on the scatterer size, with *λ* being the wavelength). The terms in the integral are all straightforward to determine: *G* is known from (2.6), *ϕ* can be taken directly from the simulation, *h* is known from the thickness map and *θ* can be calculated from geometry.

Once *ϕ*_{s} and *P*_{mn} are known, *C*_{mn} can be extracted through the use of a linear solver to invert (2.12). Clearly there are actually many scattering scenarios, giving multiple datasets for *P*_{mn} and *ϕ*_{s}; even for just a single realization of *h*(**x**), both the scattering direction and the illumination direction will vary, giving many scattering values. Also, it is desirable to have a solution which is valid across many different defects, so several defect sizes and depths will be simulated. Ultimately, (2.12) will represent a highly overdetermined system of equations, and from this *C*_{mn} can be calculated in the least-squares sense.

#### (i) Acoustic scattering with the generalized scattering model

To demonstrate the general scattering model, the acoustic case is returned to. This is simpler than a full-guided wave case, and an analytical solution exists for comparison. An elemental scatterer under the acoustic model has omnidirectional scattering, so, as discussed, only a single-order cosine expansion is needed, i.e. *M*=0. Therefore, a set of coefficients *C*_{0n} are sought to solve
*h* term retains its meaning. To test this, a frequency domain finite difference simulation is performed for several defects. This is operated at 50 kHz using parameters from the *A*_{0} mode on a 10 mm thick steel plate; node spacing was 1 mm and 240 elements were used in a circular array around the defect. All defects were axisymmetric Hann defects as described in [18], and had widths of 20, 40, 60, 80 and 120 mm, and depths of 10, 30, 40, 50, 55, 60% of wall thickness, making 26 simulations in total, including a run to obtain the incident field *ϕ*_{0}. The coefficients are then extracted from these scattered fields using the approach described previously.

There are likely to be two sources of error in the scattering behaviour predicted by these coefficients. Firstly, there is uncertainty in the modelling of the scattered fields and subsequent fitting of the coefficients to these, and secondly there is likely to be some truncation error due to fitting a polynomial to the function. It is possible to separate these and just study the truncation error. As the scattering behaviour is well known analytically, the coefficients can be determined independently, without using the scattered field
*C*_{0n} involves calculating *O*(*h*) then fitting the polynomial terms to it. The process for determining *O*(*h*) is in two steps: firstly, the phase velocity *c*, associated with the thickness, is derived by interpolating the dispersion curves, typically calculated via a tool such as DISPERSE [29], then *O* is calculated from that as in (2.4). Given the lack of simple analytical equation for the dispersion characteristics, the chosen approach is to numerically evaluate the function *O*(*h*) across a range of thickness values from 0 to 60% wall loss; the polynomial coefficients are then used to fit this function.

Figure 3 shows the comparison between the true *O*(*h*) term and two *O*(*h*) is calculated semi-analytically as discussed above, and the ‘direct’ curve compares that with its polynomial approximation calculated directly as in (2.17). Comparing these two lines first, it is clear that the polynomial approximation does introduce small differences; however, these appear small and would be unlikely to severely affect the accuracy of a reconstruction. This suggests that three polynomial terms are sufficient for this problem. The final line is calculated using coefficients extracted from the scattered field. Here, there is a slightly larger deviation at deeper depths (beyond 50% wall loss). One aspect to note is that this line and the true line match better for shallower defects, so there is a possibility that better fitting has been performed for the shallower depths. It is noted that because each defect varies from its maximum depth to zero contrast, the average depth throughout each defect will be significantly less than the maximum. Therefore, across all of the scattered fields, the majority of the scattering will come from the shallower defect values; this is why bias was given to deeper defects in the simulations. However, there are still more shallower areas in total and this could explain why better fitting is obvious in this area. Additionally, as defects get deeper, more nonlinearity tends to occur and the underlying assumptions (e.g. the lack of multiple scattering) will break down, which could cause larger errors. However, based on this study, it is expected that the scattering model with three polynomial terms should be a good approximation for the majority of defects across the range of interest.

#### (ii) Multi-modal generalization

In many scenarios, scattering occurs in a medium which can support multiple modes, and in the majority of cases some mode conversion will occur. It is noted that the current formulation could, in theory, be used to capture such multi-modal behaviour
*q* represents the incident modal wavefield with *Q* being the number of modes present. 1≤*r*≤*Q* represents the measured mode. Clearly, this formulation increases the number of parameters which need to be determined by a factor of *Q*^{2}, compared with the single-mode case. While this description is included for completeness, it has not been tested in practice and a suitable inversion approach based on it has not been derived. A very similar approach to this has been discussed in [30], and, as highlighted there and investigated in [31], a particular problem is that the K-space is not completely sampled by the combined two modes; a section of wavenumbers close to zero cannot be extracted. Another practical challenge is the complexity of experimentally measuring and exciting several different modes while maintaining the ability to separate them.

### (c) Scattering behaviour

Figure 4 illustrates the scattering patterns associated with different depth ‘elemental’ defects, given by the equation
*M*=3 and *N*=3, and are calculated by fitting the responses from defects of a variety of widths and depths down to 60% wall loss following the approach in §2b. It is observed that these results match those predicted analytically for shallow defects in [21] under Mindlin theory, which are further analysed in [30]. It is clear that in through-transmission, the two curves generally match well, with increasing deviation with depth, matching the trend seen in §2b(i), figure 3. This is more clearly highlighted in figure 4*g*, showing the transmission component as a function of wall loss, with a very similar deviation to that visible before.

In reflection, the guided wave scattering does not match the acoustic scattering behaviour well; at 10% it is larger compared with the acoustic case, but it decreases in amplitude and between 30% and 40% wall loss becomes negative. Scattering at other angles is also very different. As illustrated in figure 4*h*, there is a significant zero around 1 radian; despite the diversity of scattering behaviour with depth, at this angle the scattering amplitude remains zero for all depths. This can be problematic because the information contained in the scattered field at this angle is unavailable.

## 3. Inversion

The primary motivation for developing this model is to exploit it in inversions to generate improved images. This section begins with an overview of the diffraction tomography-based approaches for acoustic problems.

### (a) Acoustic inversion

#### (i) Diffraction tomography

Beginning from the Lippmann–Schwinger equation of (2.7)
*ϕ*=*ϕ*_{0} within the integral, and secondly that the transducer array is in the far field, so that plane wave assumptions hold; this leads to [15,32]
*A* is the far-field scattering amplitude. It is clear that this represents a 2D Fourier transform
**K**|≤2*k*_{0} can be determined. One approach is to acquire these components, interpolate them to a regular grid, then perform a 2D inverse Fourier transform to obtain *O*, as in
*K*_{x} and *K*_{y} are the Cartesian components of the vector **K**. As this provides access to components up to 2*k*_{0} this corresponds to a resolution limit of *λ*/2 [32]. However, it is possible to avoid this interpolation. If the incident and scattered vectors are expressed as functions of angles, i.e. *W* is a weighting function, corresponding to the determinant of the Jacobian matrix, to account for the change in coordinates; this is Devaney’s filtered backprojection algorithm [33]. It can be shown [34] that this weighting function has the form
*et al.* have conducted a thorough study of the behaviour of these different algorithms [36], and as explored in [37], this approach amplifies the low- and high-spatial frequency components in the image; this effect can be removed by applying a filter to generate the standard diffraction tomography image. One benefit of this approach is that beamforming images do not rely on the far-field (plane wave) assumption; as the same image can be generated regardless of transducer location, the same diffraction tomography image will be generated from near- or far-field data. As highlighted previously [8], there is often also no need for incident field subtraction with this approach, a traditional difficulty with diffraction tomography. This approach is taken as the basis for the solutions developed in this paper.

#### (ii) HARBUT

As highlighted in [14], the capabilities of diffraction tomography are limited for guided wave tomography because the typical defect sizes generate sufficient phase shift through them that Born approximation is violated (i.e. *ϕ*_{0} is not a good approximation of *ϕ*). Instead, HARBUT (Hybrid Algorithm for Robust Breast Ultrasound Tomography) [17,38] is employed. This splits the object function into two components
*O*_{b} corresponds to an approximate, known background estimate of *O*, typically estimated from bent ray tomography or a previous HARBUT iteration, and *O*_{δ} corresponds to the small, low contrast residual. As *O*_{δ} is low contrast, the Born approximation should be suitable for reconstructing it against the known background *O*_{b}. By reformulating the equations of §2a for this (see [17] for details), the scattering equation becomes
*b* is used to indicate the values of the scattered field and Green’s functions in the background, and it has been recognized that the background field can be calculated using Green’s function, assuming a point source excitation. The inversion is then performed in a similar way to before
*O*_{δ} value. Clearly the terms *G*_{b} need to be calculated; these are derived as being the same as the standard Green’s functions (2.6) except with a phase distortion accounting for the background, calculated using a fast eikonal equation solver [39]. It was shown that iterating HARBUT [10] can enable improved reconstructions to be achieved as *O*_{δ} reduces in amplitude and *ϕ*_{b} becomes a better approximation of *ϕ*. HARBUT has been shown to produce the best results based on the acoustic assumption across a range of defects [14], and will, therefore, be developed to account for the full-guided wave scattering behaviour.

### (b) Elastic guided wave inversion

Several techniques are possible, based on these methods, to derive a solution to invert guided wave data. Two approaches are explored initially, which will introduce concepts which are ultimately combined into a single algorithm to address the problem.

#### (i) Fourier component weighting

One simplifying assumption in the scattering model which could be made is that
*et al.* [21]. This solution gives *O*, and then *h* must be determined by inverting *O*(*h*) for a set of sufficiently sampled values of *h* in the range 0≤*h*≤1, which can then be interpolated for each value of *O* to give *h*.

One challenge of solving equation (3.15) is the division by **K**| values, thus enabling information which is unavailable at one frequency to be replaced with information from another. It is possible to incorporate the *O*_{ac} to a guided wave version *O*_{gw}.

This can be generalized to HARBUT. For the first iteration, the background object function *O*_{b} contains only low-frequency components, generated by a ray-tomography algorithm; these correspond to the transmission subset and as shown in §2c, the acoustic model matches the guided wave model well in this area, i.e.
*O*_{δ}, the weighting will need to be applied and this can either be applied immediately
*O*_{δ} can be filtered either immediately or at the end when they have been summed up. In the latter case, because the corrections are not applied before the corrected steering functions are calculated, there will be slight distortions to the wavefields *G*_{b}; however, it is generally the low spatial frequency components which are the most important in these calculations (and in fact the background is generally smoothed prior to calculating *G*_{b} [10]) so any differences will be negligible.

#### (ii) Iterative solution

As discussed, the assumption that the scattering behaviour is constant with defect depth is not valid for the ranges considered here. The complexity, therefore, increases significantly, because the algorithm must account for the different scattering behaviour with depth without knowledge of what the depth is. One of the most straightforward approaches would be to use a least-squares minimization based on the scattered field, minimizing a cost function defined as
*ϕ*_{s,meas} is the measured scattered field and *ϕ*_{s,model} corresponds to the modelled scattered field, as generated from (2.10), and *r*=*ϕ*_{s,model}−*ϕ*_{s,meas} is the residual. A gradient approach can be used to achieve this, and an efficient gradient calculation method can produce this derivative relative to the values of *h* in the image. This is expressed as
*ϕ*, these terms are all known. *ϕ* can be determined from any of a number of forward solutions (e.g. finite difference or FE) based on the current thickness map; probably, the fastest approach would be to follow the HARBUT formulation [17] with an eikonal solver to distort the incident field. Having calculated this gradient, a solution technique, such as the conjugate gradient method could be used to minimize the cost function by a stepping approach. While this approach has been included for completeness, it is not taken further because it requires incident field subtraction to produce *ϕ*_{s,model} and is potentially slow because the sums of (3.22) need to be calculated at every location **x** in the image for every iteration. However, the concept of least-squares fitting explored in this section can be combined with the Fourier component weighting approach discussed above, as will be explored in the next section.

#### (iii) HARBUT with an improved scattering model (HARBUTISM)

As discussed above for the separable variables case, the standard acoustic reconstruction will contain all the components of the true reconstruction, just distorted. When the variables are not separable, this still holds, but there is now no simple procedure to correct for the complex scattering behaviour as the depths are unknown. The least-squares fitting approach is again proposed, this time minimizing the cost function
**K** value in terms of plane wave illuminations and far-field measurements applied to (2.10)
*θ* will be a function of **K**.

The gradient is calculated, as before, relative to *h*(**x**)

This approach will be referred to as HARBUTISM, HARBUT with an Improved Scattering Model. There are a number of advantages to this approach. It is much faster; the 2D Fourier transform significantly reduces the time taken compared with the multiple sums required previously. The approach also means that the advantages of HARBUT are maintained: it can generate accurate images with a wide range of defects and does not require the incident field to be subtracted. Also the singularity problem discussed for the separable variables case is minimized because division by the scattering pattern is not required. Instead, the two fitted functions will both have low-amplitude components at particular values of |**K**| corresponding to the angle where scattering reaches zero; however, the issues arising from this can be compensated for via regularization approaches or combining multiple frequencies as mentioned previously.

In this paper, regularization is applied at each iteration of the gradient descent algorithm to aid convergence. Two forms of regularization are used: positivity and filtering. Positivity exploits the prior information that it is not possible to have negative wall loss, i.e. the plate thickness cannot increase, only decrease. This provides a restriction on what the wall loss values can be; therefore, at each iteration, the image is forced to match this restriction by setting any values *h*<0 to 0. Filtering can remove spurious unwanted high spatial frequency components. Filtering is applied within HARBUT anyway, as part of the beamforming to diffraction tomography filter [37]; this removes any wavenumber components above the fundamental 2*k*_{0} limit. In fact, the filter cut-off values are selected to be slightly lower than this, to 1.7*k*_{0}, sacrificing resolution to improve stability. This filter is also applied in HARBUTISM at every iteration in order to minimize high spatial frequency artefacts.

It is noted that gradient descent methods as proposed here can suffer from convergence to local, rather than global, minima. However, the starting point of the iterative process is the HARBUT reconstruction, which has been shown to generally be close to the true solution. Analysis in this paper will evaluate several simple cases and two complex defects, providing a degree of confidence that local minima do not cause a problem; however, future work may evaluate the stability of the algorithm across a much larger set of scenarios.

## 4. Results

The test cases presented here are similar to those used previously in [14,18], so only the key parameters are given here. In all examples, the *A*_{0} mode has been used at 50 kHz on a 10 mm thick steel plate, although results are normalized to be valid for other thicknesses, provided the frequency is adjusted to ensure the same frequency–thickness product.

### (a) Simple defects

The configuration is illustrated in figure 5. As shown in figure 5*a* the array surrounds the defect and has a radius of 180 mm, with 240 transducers. Axisymmetrical Hann-shaped defects were used, with thickness functions defined as
*D* and width *W*. Owing to symmetry, only a single illumination was needed, and the information from that could be replicated to make a complete dataset containing all the send–receive combinations. The FE software package Pogo [28] was used to simulate the different defects, with a uniform structured mesh squashed to model the defects as shown in figure 5*c*. At the frequencies used, only the *A*_{0} and *S*_{0} modes are present, so the antisymmetric source shown in figure 5*d* was used to excite pure *A*_{0}. Data from these simulations are included as electronic supplementary material.

Figure 6 presents the results for defects of different widths and depths. Overall, it is clear that the accuracy is improved using the improved scattering model in HARBUTISM versus relying on the acoustic to velocity mapping approach of HARBUT. This is particularly notable for narrower defects, indicating a resolution improvement; this fits with the observed resolution limit of the velocity-to-thickness mapping approach of around 2*λ* first identified in [14]. With the deeper defects, larger errors are typically visible. The polynomial approximation will inevitably contain more errors for the larger values of *h* where the truncation errors become larger; this will introduce errors both when calculating the coefficients *C*_{nm} and when using them in the inversion scheme. HARBUTISM can generally be expected to perform best for shallower defects.

#### (i) Experimental comparison

Data from the two axisymmetric experimental cases of [14] were used for experimental validation. The defects had widths of 120 mm and 60 mm, and depths of 6 mm and 3 mm, respectively. Measurements were taken around half the array with a laser vibrometer at 129 locations, as illustrated in figure 5*a*, and replicated to give a full array, ultimately forming an array of 256×256 measurements.

Figure 7 presents the experimental results, with the equivalent simulated versions repeated for comparison. For both defects, it is clear that the experimental data show the same trends visible in the simulations. For the smaller defect, HARBUTISM clearly avoids the resolution limit present in HARBUT. In both experimental datasets, HARBUTISM tends to overshoot more than the simulated equivalent, but this is generally by a small amount. It is likely that this has come from experimental uncertainties, such as potential anisotropy, present in the plate. Another potential source of errors is the data replication process. Because of this, noise, which would normally be incoherent between measurements, and hence be largely cancelled in the final image, is now coherent and forms axisymmetric ring artefacts as seen in [10]. It could be expected that such artefacts could also be superposed onto the defect reconstruction itself, affecting the accuracy.

### (b) Complex defect

The initial tests have demonstrated the performance of HARBUTISM for both simulated and experimental defects. However, it is important to test its behaviour when being extended to a fully complex defect. For these purposes, a laser scan of a real corrosion defect is used. This is scaled to make its maximum depth 50% of wall thickness, and the resulting thickness map is illustrated in figure 8*a*. As before, a full FE simulation is performed; this time, however, it is necessary to simulate all 240 sources in order to fully capture the field as there is no symmetry.

The reconstruction using the standard velocity–thickness mapping, calculated with HARBUT, is shown in figure 8*b*. Clearly significant detail is missing from this reconstruction compared with the true thickness map. When full-guided wave scattering is accounted for with HARBUTISM, figure 8*c*, more of the finer features are visible in the reconstruction. Importantly, the deepest part of the defect is now a lot closer to being accurately captured by the imaging method. This is validated in the cross section through the deepest point of the defect, shown in figure 8*d*. From this, it is clear that the error between HARBUTISM and the true depth is less than 10%, while the HARBUT error is closer to 18%. This is largely caused by the differences in resolution between the two methods.

Resolution can be analysed by looking at the spatial frequency components. These are obtained by performing a 2D discrete Fourier transform on the values of *h* for the true, HARBUT and HARBUTISM images of figure 8; the results are shown in figure 9. It is clear that HARBUTISM extracts significantly more spatial frequency components than standard HARBUT. Quantifying resolution is challenging from such plots, but it is estimated that HARBUT can extract components in the range |**K**|<0.6, while HARBUTISM achieves around double this, |**K**|<1.2, corresponding to a resolution limit of less than *λ*. This is worse than the theoretical limit |**K**|<2; in part this is caused by the regularization filter removing components above |**K**|=1.7, but also likely to be related to the challenges of accurately extracting and imaging from the very weak components scattered in the reflection subset of the data. Future work may examine approaches to reach this limit.

A further example has been tested to examine this. This is shown in figure 10*a* and consists of several axisymmetric Hann shapes superposed. As before, the HARBUT resolution is limited (figure 10*b*), while HARBUTISM is more accurate (figure 10*c*). The spatial frequency components of figure 10*d*–*f* shows similar results to those obtained before, with HARBUT accurately capturing details up to around |**K**|<0.6. However, for this defect, HARBUTISM appears to produce a higher resolution reconstruction, up to around |**K**|<1.6, close to the regularization filter level. It is unclear at present exactly why this case achieves better results. Notable differences include the shape being simpler than the real corrosion patch, the in-plane extent being smaller, and also the depth being slightly less than before, and it may be one of these which is related to the improvement.

One feature to note in figure 10*e* is the darker circle marked with the arrow, which is not visible in either of the other spatial frequency spectra. This occurs at a constant value of |**K**| and, therefore, corresponds to a particular scattering angle. This angle is an angle at which no scattering occurs from the guided waves, and matches the zero identified in figure 4*h*.

#### (i) Speed

HARBUTISM has not been highly optimized for speed; potentially there could be good savings from using the high parallelism of GPU (Graphical Processing Unit) technology for this task. However, even in its current form, written in Matlab, it runs at an acceptable speed. For the initial complex case, running on an HP Z820 workstation with four 8-core Intel E5-2665 Xeon processors operating at 2.4 GHz, the ray tomography time was 2.3 s, standard HARBUT took 12.7 s in addition to this, then HARBUTISM took a further 18.0 s, making a total of 33 s to generate a full image of a complex defect. This indicates that the system can deliver accurate images in a practical time-frame.

An alternative approach to solving this problem is to use a full numerical model to generate data, then iteratively update it until the numerical data match the measured data. This could be done via a gradient descent approach or with one of the many approaches to sample the parameter space in order to find a global minimum. Either approach requires the full model to be run many times to match the two datasets. However, a complete forward model, even on multiple GPUs with Pogo [28], typically takes many hours or days to produce a full dataset, and this must be repeated many times to achieve convergence. Even with the rapid pace of computational improvements, it is unlikely that this approach will be a practical solution for guided wave tomography in the foreseeable future.

## 5. Discussion

### (a) Density effect

A recent paper [18] discussed the effects of density in guided wave scattering, and highlighted how it can cause critical errors in the reconstruction if neglected. It is noted that this is the effective density that the guided wave experiences, rather than the material density; this effective density term is proportional to thickness. The improved scattering model has so far not explicitly accounted for density, so the effects of it will now be explored.

The object function incorporating density is defined as [40]
*ρ*(**x**) is the local density. For simplicity in this analysis, it is assumed that thickness (and hence density) changes are small and that *ρ*^{1/2}(**x**)=*ρ*_{r}(**x**)=*ρ*_{r0}+*δρ*_{r}(**x**), where *ρ*_{r} has been introduced to signify the square root of density. The density component of the object function then becomes
*ρ*_{r0}+*δρ*_{r})^{−1}. As, under the far-field Born approximation, the scattered field can be mapped to the Fourier transform of the object function, the density component of scattering can be written as
*θ* as
*δρ*_{r} term can be approximated easily by a polynomial in *h*^{n}. This behaviour can, therefore, be captured by the generalized scattering model, and will automatically be accounted for in the coefficients extracted from the simulation. This integrated technique is arguably a more elegant solution for density correction compared with the previous approach [18] where an initial thickness map was generated ignoring density, then iteratively updated with a correction term based on the implied density.

This analysis has been performed using the assumption that the thickness variations were small, allowing the equations to effectively be linearized. However, for the deeper defects, it is likely that the density scattering behaviour will be similar, and, therefore, the improved scattering model should still be able to capture the behaviour, provided it has enough coefficients. It is recognized, however, that imperfect fitting of this could be a cause of additional errors seen for the deepest defects.

### (b) Stability

Any algorithm which extracts more information from the same data is likely to suffer more when the data have errors in it. The use of experimental data has helped to provide validation of the performance of HARBUTISM in non-ideal scenarios, dealing with both noise and systematic errors; however, it is challenging to quantify the effect of such errors. The performance with different errors will be evaluated in a more controlled manner here using synthetic data.

Systematic errors are considered first. These errors are coherent between different signals causing them all to be distorted in the same way. Examples of this could be consistently mispositioning transducer locations or a change in temperature causing a change in the sound speed. To study this, the transducers will be effectively offset; however, to avoid re-running the full FE simulation, the modelled transducer locations are kept the same, but their locations are mis-reported to the imaging algorithm, which will have an equivalent effect. The reported transducer locations are scaled by a constant factor in both in-plane dimensions, making the circular array slightly bigger or smaller. It should be noted that this is a slightly artificial scenario as, in general, the locations will be calibrated beforehand (e.g. [41]) to ensure the effect of any such errors are minimized in practice. Figure 11*a*–*e* shows the results from this scaling in a range of ±2%; for reference, a 2% error corresponds to around 20% of a wavelength for the through-transmission ray path, or 1.3 radians phase shift. It is clear that in this range the shape still remains well defined, and HARBUTISM, while inevitably being sensitive to the errors to an extent, does remain stable. When the given transducer array size is reduced, the algorithm reduces the reconstructed sound speeds to compensate, which corresponds to the background errors clearly visible in figure 11*a*, and less so in figure 11*b*. When the size is increased, at 1%, a slight underestimate in wall loss is apparent (corresponding to a tendency to overestimate sound speed). At 2% the HARBUTISM fitting routine has not worked because of the errors in the data preventing a better solution from being found, and the standard HARBUT reconstruction, which is the starting point of the fitting routine, is produced.

In many scenarios, the effect of random noise may be significant. This can be caused by electrical noise in the acquisition system, or by any other incoherent transient effect on the system. To test these effects, noise is added to the measured data matrix prior to imaging. This is done in the frequency domain at the frequency of interest, and has normally distributed amplitude, with a standard deviation at a particular fraction of the true matrix in order to give a certain signal-to-noise ratio (SNR). The phase of each component was uniformly distributed from 0 to 2*π*. The results are shown in figure 11*f*–*i* for 6 dB, 10 dB, 15 dB and 20 dB, respectively. It is clear that as the SNR increases, so does the image quality. HARBUTISM does show quite significant, large artefacts for the low SNR values, but by the 20 dB point (which most modern acquisition systems can achieve easily) the image shows very little difference from the noise-free version.

### (c) Sampling

A full analysis of sampling is considered beyond the scope of this paper; however, it will be briefly discussed here. Provided the wavefield is fully sampled (i.e. with less than *λ*/2 transducer spacing) the Born-based inversion will have sufficient data. For the array radius of 180 mm and the wavelength of approximately 40 mm, the minimum number of transducers, therefore, needed is around 57; going beyond this with over 200 transducers, as has been done in this paper, will have no effect on the resulting image. If the scattered field is used to image (rather than the total field used here), then the sampling can often be reduced further [42]. Future work will study the case when fewer transducers are used and subsampling occurs.

## 6. Conclusion

The majority of guided wave tomography approaches use a velocity-to-thickness mapping approach which assumes that elastic guided waves scatter from corrosion defects in the same way as acoustic waves scattering from velocity variations. This results in important inaccuracies in the reconstructions. This paper has introduced HARBUTISM, a technique which avoids this assumption, more accurately capturing the guided wave scattering when producing the thickness map. This leads to significantly improved accuracy in the reconstruction, which is partly visible as an improvement in the resolution of the image. From a resolution limit of standard velocity inversion of around 2*λ*, HARBUTISM has been shown to achieve at least *λ* and is theoretically capable of *λ*/2. Such resolution improvements are critical for fully capturing the complex peaks often present in real corrosion defects, which is necessary to enable the maximum wall loss to be accurately extracted. This is very important to enable the petrochemical industry assess the life of their components.

HARBUTISM was tested and compared with the velocity-to-thickness approach for a number of cases: simple axisymmetric shapes with simulated and experimental data, and on a complex simulated example taken from a laser scan of a true corrosion defect. Accuracy was improved in the majority of cases. However, for some deeper defects where wall loss was beyond 50%, modest overestimates in wall loss of around 5–10% of wall thickness were visible, which is generally not significant. It is also noted that by overestimating damage, this provides a conservative estimate of remaining life. HARBUTISM was tested with artificial sources of errors and noise, confirming itself to be robust, and produced these images in practical time-frames of less than one minute.

Finally, it should be noted that the improved scattering model is a general approach. Provided that scattering data exist to which the coefficients can be matched, it is then possible to perform an inversion to produce images. It is, therefore, expected that this approach, in combination with modern high-speed simulation tools, could be applied to other areas where the inversion of complex scattering behaviour is needed.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Competing interests

I declare I have no competing interests.

## Funding

This work was funded by EPSRC under grant no. EP/M020207/1.

## Acknowledgements

The author thanks Prof. Peter Cawley for reading the manuscript.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3517530.

- Received August 19, 2016.
- Accepted October 5, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.