In corrosion assessment, ultrasonic wall-thickness measurements are often presented in the form of a color map. However, this gives little quantitative information on the distribution of the thickness measurements. The collected data can be used to form an empirical cumulative distribution function (ECDF), which provides information on the fraction of the surface with less than a certain thickness. It has been speculated that the ECDF could be used to draw conclusions about larger areas, from inspection data of smaller sub-sections. A detailed understanding of the errors introduced by such an approach is required to be confident in its predictions. There are two major sources of error: the actual thickness variation due to the morphology of the surface and the interaction of the signal processing algorithm with the recorded ultrasonic signals. Parallel experimental and computational studies were performed using three surfaces, generated with Gaussian height distributions. The surfaces were machined onto mild steel plates and ultrasonic C-scans were performed, while the distributed point source method was used to perform equivalent simulations. ECDFs corresponding to each of these surfaces (for both the experimental and computational data) are presented and their variation with changing surface roughness and different timing algorithms is discussed.

## I. INTRODUCTION

Corrosion is a significant problem in the oil and gas industry. Estimates put its cost to the petroleum sector at around $8 billion annually in the USA alone.^{1} To track the progress of corrosion, the regular inspection of vulnerable infrastructure is required. Inadequate monitoring of the progress of corrosion can lead to catastrophic failure of the plant^{2} along with severe safety, environmental, and economic consequences, and potential criminal prosecution.^{3} Inspections are performed using non-destructive testing (NDT) techniques, at specified time-intervals. A brief overview of these techniques can be found in Ref. 4 with ultrasonic thickness measurement being the most commonly used.

Ultrasonic tests are often performed as a C-scan. C-scans are used to construct area color maps of the measured thickness. At every point across the inspection area, a thickness measurement is taken. An image of the measured thickness distribution is formed by representing each measurement with a colored patch, the color or gray scale chosen to be representative of the measured thickness.

An alternative method of presenting C-scan data is a cumulative wall thickness distribution function (CDF). Many examples of CDFs, extracted from C-scans obtained from pressure equipment with in service degradation, are presented by Stone.^{5} These wall thickness distributions display regular and ordered behavior, showing an exponential tail and an overlying Gaussian profile. The exponential tail is associated with more localized corrosion (very few, deep defects), while the Gaussian profile is attributed to general corrosion and/or the as-built thickness variation across the entire surface.

It has been suggested that CDFs could be used for partial coverage inspection (PCI) as part of non-intrusive inspection approaches;^{6} the use of CDFS based on corrosion mapping data has been included in an industry recommended practice.^{7} An inspector takes measurements across an accessible area, which is assumed to be under the same conditions as an inaccessible or uninspected region. The CDF calculated from the inspection can then be used to make an assessment of the remaining area. In order for this to be effective, an understanding of the errors associated with C-scans of rough surfaces is required.

Ultrasonically measured thickness measurements are a combination of the interaction of the ultrasonic pulse with the corroded surface and any noise introduced by the signal processing algorithm and the signal acquisition system. Jarvis and Cegla investigated the stability of three commonly used timing algorithms used to extract wall thickness measurements from signals collected using a permanently installed shear wave monitoring system.^{8} It was found that the wall thickness estimate changes significantly for different instances of a rough surface (with the same statistical description), if a different algorithm is used. The question has been raised as to whether the longitudinal waves used by standard transducers will be affected in the same way. It is the aim of this paper to describe the effect of the surface roughness and the choice of timing algorithm on ultrasonic thickness measurements taken using a longitudinal probe.

We begin by describing a typical C-scan set up, methods of generating known rough surfaces and commonly used timing algorithms (Sec. II). A detailed analysis of the model used for the ultrasonic simulation follows (Sec. III), with comparison to an analytical solution, convergence studies, and a study of boundary conditions. Experimental validation of the model is then described (Sec. IV) and the results are presented in Sec. V.

## II. BACKGROUND

### A. Corrosion mapping

The corrosion mapping set-up consists of a 6 mm diameter 5 MHz longitudinal wave transducer coupled to the surface using water (or another suitable couplant). There are two alternatives to couple the transducer to the part: directly placing the transducer on the surface of the part with a small amount of couplant (a contact scan) or placing the part in a water bath. In industry most scans are performed as contact scans as it is infeasible to use immersion scans for in-service pipework. For the purposes of this paper, the simulations model a contact scan as it is less computationally expensive to model than an immersion scan. In contrast, the experiments use an immersion scanning set-up to ensure consistent coupling between the transducer and the plate across the surface of the part.

For the simulations, the front surface is assumed to be in good condition, so that the transducer makes a flat contact (Fig. 1). To make a thickness measurement, an ultrasonic beam is then radiated into the steel wall, which is reflected from the internal, corroded, rough surface and recorded (Fig. 2). The thickness is calculated from the time of flight of this pulse. In real situations, multiple reflected pulses will be received. However, the signal processing is restricted to the first reflection for the purposes of our study, as it is often used for thickness measurements.

To obtain a map of wall thickness measurements, the transducer is moved a small distance in either the *X* or the *Z* directions, collecting an ultrasonic signal at each point across the inspection area. The beam profile at the backwall is collimated, so each measurement will only probe a small area directly under the transducer. The reflected ultrasonic signal from this patch will consist of scattered energy from the surface roughness (Fig. 2). Consequently, the signal shape can change substantially between closely separated measurements. These two effects lead to a natural variation in the thickness measurements across the surface, determined by the characteristics of the transducer and the surface roughness.

### B. Surface roughness

The surfaces studied in this paper are generated by a large number of independent, random events. Although the height probability distributions of these events may be very different, the overall distribution will tend to a Gaussian.^{9} There is strong evidence in the literature to support this claim, for surfaces generated by general uniform corrosion.^{5,9–11} For more localized corrosion, the height distribution can follow exponential distributions.^{5,12,13} The type of corrosion that can be expected depends on the damage mechanism that is most likely to occur in the vessel, for the internal operating conditions (temperature, pressure) and chemical conditions (pH, species present). For the purposes of this paper, we restrict ourselves to surfaces generated by a general uniform corrosion mechanism.

To demonstrate that a Gaussian profile is representative of the type of damage which should be expected in the case of general corrosion, surface profile measurements were taken from a pipe sample. The sample was retired from a unit exposed to high temperature sulfidation corrosion (a uniform corrosion mechanism). A Talysurf^{TM} surface profile measurement instrument^{14} was used to measure the profile. Talysurf^{TM} draws a stylus across the surface of the material, measuring the deflection of the stylus with an interferometer, extracting the surface profile. Several height profiles from around the inner radius of the pipe were measured. An example of a histogram of the measured heights is given in Fig. 3.

Simulation of multiple rough surfaces was performed using the algorithm described by Oglivy.^{10} Other approaches exist,^{15} however, Oglivy's was chosen due to it is computational simplicity. Oglivy's method generates a set of uncorrelated random numbers and performs a moving average, producing a set of Gaussian correlated random numbers. The probability distribution function of such a surface is given by

where $h$ is the height of the point above the plane of the surface and $\sigma $ is the root-mean-square (rms) surface variation, which controls the vertical extent of the roughness.

The local heights are correlated, in the horizontal direction, with a Gaussian weighted correlation function

with a correlation length $\lambda c$, which controls the distance at which the heights of two points becomes statistically independent. The script used to generate the rough surfaces is the same as in Ref. 8 and implemented in matlab.^{16} Full details of the roughness generation algorithm can be found in Ref. 10.

Surfaces which have different correlation lengths in different directions do occur in practice. However, we restrict ourselves to an isotropic correlation length. Methods to generate non-isotropic correlation profiles can be found in Ref. 15.

### C. Timing algorithms

The transducer operates as a pulse-echo sensor, with the only scattering occurring at the pipe boundaries. The thickness of the component is calculated as

where $t2$ is the time of arrival of the pulse reflected from the rough surface, $t1$ is the time of arrival of the reflection from the front of the component, and $vL$ is the speed of sound. The time of flight (TOF) is defined as $t2\u2212t1$.

There are many ways of measuring the time of flight. The performance of three common timing algorithms, presented below, is compared in this paper. In all of the algorithms, the signals were interpolated, to increase the accuracy of the TOF measurement.

#### 1. Envelope peak detection

The Hilbert transform was used to calculate the envelope of the pulse. The time at which the maximum occurs in this envelope is $t2$,

where $f(x)$ is the reflected pulse, $x$ is an integration variable, and $t$ is time. Similarly, $t1$ is the time at which the maximum of the envelope of the transmitted pulse occurs.

#### 2. Cross-correlation

The outgoing pulse *g*(*t*) and the received pulse *f*(*t*) are cross-correlated. For $J$ samples of the signal, the cross-correlation at sample $k$ is given by

The minus sign in front of $f(t)$ accounts for phase reversal on reflection. The maximum point in $s(k)$ is the time of flight,

where $t(max(s(k)))$ is the time at which the maximum value in $s(k)$ occurs. In the experimental part of the study, the outgoing pulse is the reflection of the ultrasound from the flat front surface.

#### 3. Threshold first arrival

The Hilbert transform was used to calculate the envelope of the pulse. The amplitude of the envelope is normalized relative to its maximum amplitude.

A threshold amplitude was selected and the time when the part of envelope corresponding to the reflected pulse, crosses the threshold is $t2$. $t1$ is the point at which the envelope of outgoing pulse crosses the threshold.

## III. MODELING REFLECTED ULTRASONIC SIGNALS

To model a C-scan, it is necessary to simulate a large number of thickness measurements. For a typical scan of $200\u2009mm2$ of a component with $dx=dz=1\u2009mm$, the number of thickness measurements taken is 40 000. Therefore, a fast simulation method is required in order to complete the modeling in a reasonable amount of time.

Finite element methods (FEM) are the standard for simulating ultrasonics. Commercially available 3D FEM codes are too computationally expensive for our purposes and recent developments, such as pogo,^{17} are restricted by the memory available to the latest graphics processing units. A 2D to 3D FEM conversion has been suggested,^{8} but this relies on a specific transducer geometry.

As we are using a zero-degree longitudinal wave transducer, it can be assumed that mode conversion has little effect on the received signal. Therefore, the distributed point source method (DPSM), which models acoustic fields, can be used to reduce the computational cost of ultrasonic modeling. It has been show that DPSM performs much faster than standard FEM codes for the types of situation considered in the current work.^{8}

### A. The DPSM

#### 1. Frequency domain calculations

The DPSM is a mesh-free semi-analytical technique which offers a fast alternative to the finite element method for acoustic field calculations. The DPSM consists of the modeling the transducer with a set of active point sources^{18} and the rough surface with a set of passive point sources.^{18} These point sources are offset from the exact position of the transducer and the rough surface by a small distance determined by the spacing between the point sources.^{19}

As the pressure Green's function of a single point source is known,^{8} the amplitudes of the point sources can be calculated to satisfy the boundary conditions at the transducer face and the rough surface.^{20} Once these amplitudes have been calculated, the pressure at any point in the simulation cell is a sum of the pressure from each point source,^{21} which can be expressed as a matrix multiplication.^{22} The DPSM has been applied to a number of problems in acoustics.^{8,23–25} Readers interested in more detail are referred to Placko *et al.*^{26}

As there is no commercially available software for DPSM modeling, a software package was developed for the purposes of this study. The program performs the calculations in parallel using mpi,^{27} handles the Fourier transforms with the fftw library,^{28} and the boost library was used to handle various other mathematical routines.^{29} For all of the simulations in this paper, the geometrical set-up of the DPSM model is shown in Fig. 1. Active source points (the transducer) are placed on the plane *Y = *0, with passive point sources (the rough surface) placed with a mean plane of $Y=25\lambda /3$.

Figure 4 shows a comparison between the DPSM solution and the analytical solution derived by Mellow^{30} for the pressure on the axis $(0,\u2009Y,\u20090)$ radiating from a 5 MHz 6 mm diameter longitudinal transducer into an acoustic medium with the same density and compressional wave velocity as steel (i.e., $vL=5960\u2009m/s\u21d2k=5.27\xd7103\u2009rad\u2009m\u22121$). The discrepancy shown in Fig. 4 has been well-documented and has been attributed to poor matching to the uniform pressure boundary condition at the transducer face due to the limited number of point sources.^{31} Alternative formulations of the DPSM have been developed with the goal of improving the matching to the uniform pressure boundary condition and they have shown some success at reducing this discrepancy.^{32,33}

If the uniformity of the boundary conditions was the cause, then increasing the number of active point sources should reduce this discrepancy. However, despite further calculations of the on-axis pressure, with an increasing active point source density, this discrepancy remains (Fig. 4). The analytical solution uses two boundary conditions: the uniform pressure boundary condition, where the pressure on the transducer face is $p(x,\u2009y,\u2009z)=1$; and the zero pressure boundary condition, where everywhere on the plane (*x,0,z*), $p(x,\u20090,\u2009z)=0$, with $x2+z2>a2$.

The zero-pressure boundary condition is not explicitly satisfied by the DPSM models. For the analytical solution, a cross-section of the pressure distribution across its surface, would show uniform pressure (equal to unity) across the transducer. At its edges, the pressure drops immediately to zero. In comparison, the cross-section of the pressure distribution calculated by DPSM, tapers off slowly (Fig. 5), leading to a “halo” around the outside of the transducer, where the zero-pressure boundary condition is not fulfilled.

To better satisfy the zero-pressure boundary condition, one can add a number of boundary source points around the outside of the transducer. This can be achieved by simply increasing the size of the propagation matrix^{34} and the boundary condition vector.^{22} For *M* active source points and *K* boundary points:

where *S′* is the set of source points and boundary points and each element of $QTS\u2032$ is the pressure Green's function between each target and source point.^{8} The boundary condition vector is given by

where $0K\xd71$ is a *K* dimension vector of zeros. The boundary and active point source amplitudes can then be calculated inside the DPSM framework.^{20}

Zero-pressure boundary points reduced the pressure around the edge of the transducer to close to zero (Fig. 5, solid black line), with the addition of four layers of boundary points (with a point separation of $\lambda /12$). While the pressure conditions are much more closely matched, there are still some differences in the pressure profile produced by the DPSM and the profile used for the analytical solution. This is a consequence of trying to model a continuous disk by a number of point sources.

The DPSM pressure profile radius in Fig. 5 (solid black line) is slightly larger than the prescribed transducer radius (3 mm), which would lead to phase differences in the near field and amplitude discrepancies in the far-field. An effective radius for this pressure distribution can be calculated by minimizing the difference between the DPSM and analytical solutions across the axis with respect to the radius. This leads to close agreement between the analytical and the DPSM solutions (Fig. 6).

#### 2. Time domain

Common practice in wall-thickness measurements is to use time-domain pulses. The DPSM calculations outlined so far are performed in the frequency domain and require extension for time domain pulse propagation. The method used here follows Jarvis and Cegla.^{8} Alternatively, the time-domain Green's function could be used, although this approach is unsuited to long pulses.^{35,36}

The frequency-domain representation of the outgoing pulse, *f*(*t*), is calculated at the transducer face,

where *F*($\omega $) is the frequency domain of the incident pulse, $\omega $ is the angular frequency, and $t$ is time. A frequency DPSM calculation is then performed for each frequency component.^{22} A single reflection is modeled and the reflected frequency domain at each point modeling the transducer face is constructed. The arithmetic mean of all of these frequency domain representations, gives the frequency domain representation of the received signal. The time-domain received signal is calculated by performing an inverse Fourier transform,

The outgoing pulse used was a 5 MHz Hanning windowed toneburst, with five cycles. To reduce the computational time taken for these simulations, any frequency component with an amplitude of less than 1% of the largest component was set to zero. These components were found to have little effect on the overall signal.

To ensure the accuracy of this approach, the analytical solution was used to calculate the pulse propagated to a single point at $(0,\u2009503\lambda ,\u20090)$ which was compared to the equivalent DPSM signal. The frequency domain representation of the outgoing signal was calculated using Eq. (9). At each frequency the amplitude at $(0,\u2009503\lambda ,\u20090)$ was found using the analytical solution and used to form the frequency domain representation of the propagated pulse. This representation was then transformed into the time domain using Eq. (10), which gives the analytical received signal at $(0,\u2009503\lambda ,\u20090)$. Close matching between the DPSM (using an active point source separation of $\lambda /12$) and analytical signals was achieved (Fig. 7), with a 2% maximum difference of the Hilbert envelopes of the signals.

### B. Validation by modeling the resilient disk

As this paper uses a relatively new method for calculating the time domain pulses, a small convergence study has been performed to select appropriate values for the passive point source density and the dimensions of a patch of rough surface.

#### 1. Passive point source density

An analytical solution for a reflected signal from a flat backwall can be calculated by reversing the phase of the analytical signal. This models the total distance traveled by the pulse and its reflection from a flat backwall $253\lambda $ away from the transducer.

The DPSM was used to calculate the reflected signals from four backwalls with different point source densities. The maximum error on each signal was found by calculating the percentage difference between the maximum of the DPSM solution and the analytical solution. The passive point source density when the error is reduced to 1% or less, was used for the calculations. This was found to be $147\lambda \u22122$, corresponding to a passive point source separation of approximately $\lambda /12$ (Fig. 8).

However, this paper is concerned with rough surfaces. Therefore, the parameters need to be checked for a surface of varying height. A surface with a sinusoidal height variation was used as a model surface. Its amplitude and wavelength were chosen to be of a similar extent to rms and correlation lengths. The reflected signals from this surface with various point source separations were calculated using DPSM.

For this case there is no available analytical solution for the reflected signals. Therefore, the maximum error was found by calculating the percentage difference from the signal calculated with the largest passive point source density. The calculation was taken to be converged when the reduction in error was less than 1% upon doubling the number of point sources per unit area. This was achieved for a point source density of $147\lambda \u22122$ $[100\u2009mm\u22122]$ (circles in Fig. 8).

#### 2. Backwall size

The active and passive point source densities found in Secs. III A and III B 1 were then used to determine the minimum patch size required. The reflected signal from square patches ranging from $5\u221211\u2009mm$ ($25\lambda /6$ to $45\lambda /6$) in dimension were calculated using the DPSM and compared to the analytical signal. The maximum error was calculated in the same way as in Sec. III B 1 and reduced to less than 2% with a square patch size of 9 mm × 9 mm (Fig. 9), so this was chosen as the patch size for the C-scan model.

### C. Modeling a C-scan

The previous sections showed which parameters were required to ensure DPSM simulations with an error of less than 2% in the maximum amplitude of the Hilbert envelope of the signal. These were used to simulate signals collected from a wall thickness C-scan. The beam from the transducer is very collimated, only probing a small footprint directly under the transducer. Therefore, each simulated measurement only needs a small region of the surface to accurately model the reflected signal. This allows for the C-scan to be split into a large number of independent measurements, which makes the problem computationally tractable.

Each measurement is represented by a single simulation cell. A cell consists of the transducer placed on the (*X*,0,*Z*) plane (as in Fig. 1), with a corresponding patch of the rough surface placed in the (*X*, $253\lambda $, *Z*) plane. The active point sources which model the transducer are surrounded by four layers of boundary point sources at which the pressure is set to zero (as in Sec. III A 1). The patches of rough surfaces were taken from three different 200 mm square surfaces.

The surfaces were generated with a Gaussian distributed height profiles, with rms surface variations of $\lambda /12$, $\lambda /6$, and $\lambda /4$, and $\lambda c=2\lambda $. These surfaces were split up into 9 mm × 9 mm ($14\lambda /2$ × $14\lambda /2$) patches, with a lateral separation of 1 mm. These patches each correspond to a thickness measurement taken 1 mm apart, in a square grid.

Ultrasonic simulations for each cell were performed, producing a recorded signal. Thicknesses were calculated from each cell and used to calculate CDFs (Sec. V B) for the plates.

## IV. EXPERIMENTAL VALIDATION OF THE MODEL

### A. Machining a roughness profile

To further validate the model, 200 mm square patches of roughness (rms = $\lambda /12$, $\lambda /6$, and $\lambda /4$, and $\lambda c=2\lambda $) were machined onto three 300 mm square mild steel plates. To manufacture the plates, roughness profiles were generated using the script in Sec. II B and used to build a Solidworks^{TM} (Dassault, Vlizy-Villacoublay) model of the plate. The height profile was machined onto the plate using a Bridgeport^{TM} CNC machine (Series 2, Interact 4, Bridgeport, New York).

The minimum radius of curvature achievable by the CNC machine was 2 mm. Any surface features with a radius of curvature less than 2 mm were filtered from the point cloud data, by performing a spatial Fourier transform of the surface and removing the relevant frequencies. Furthermore the rough surfaces taper down gradually to a flat surface, for the last 2 mm of the rough patch, to allow cutter access.

To validate the surface variation of the plates, the Talysurf^{TM} (Taylor Hobson Ltd, UK) surface profilometer was used to measure the roughness profile of several lines across the surface of each of the plates. From these profiles, the rms and correlation length of the surface was calculated. Figure 10 shows the measured rms surface variation against the target rms surface variation. There is good agreement between the rms of the real surfaces and the rms of the computer generated surface. Furthermore, Fig. 10 shows that the measured correlation lengths are all within 0.2 mm of the target correlation length.

### B. Experimental set-up

The experiments were performed using the set up in Fig. 11(a). An Ultrasonic Sciences Limited USL Scanner, with an Olympus NDT 5 MHz 6 mm diameter plane longitudinal transducer (near field distance 10.8 mm), was used to scan the plates containing the rough surface. The scanning frame arm draws the transducer across the surface of the plate, sending an ultrasonic pulse at pre-set positions (a 1 mm × 1 mm resolution raster scan).

The pulse travels through the water, entering the plate at a distance of 32 mm from the transducer face. Then, the pulse travels through the metal, until it reflects off the rough surface. The scattered signal is then recorded with the same transducer in the form of a time-trace (A-scan). The thickness corresponding to each time-trace was extracted in the post-processing using the previously described timing algorithms.

The rough surface is in the far-field of the transducer. As a consequence the field at the backwall will be of a very similar shape to that in the simulations, albeit with a reduced amplitude.^{37}

## V. RESULTS

### A. Signals

The time domain signals collected from across the surface give insight into the measurement process. The transmitted pulse travels through the steel plate and is reflected from the rough surface. Figure 11(b) shows a raw A-scan signal acquired by the equipment. It shows the front reflection from the flat surface of the steel plate and the reflection from the rough surface. We were only interested in the reflection from the backwall and therefore the time traveled through the plate was computed by setting the time at which the pulse enters the plate as time zero. The distance traveled in the steel plate was calculated by multiplying the arrival time of the reflection by the longitudinal wave velocity in steel (5960 m/s). Subsequently traces containing only the backwall reflection will be shown (see Fig. 13)

The signals in Fig. 12(A) were collected from the rms = 0.3 mm plate (using the immersion scanner), from points on the surface a small distance apart (1 mm). There is a clear difference between the two signals. For the solid signal, there is constructive interference between components of the signal reflected from different parts of the rough surface, leading to a large signal amplitude. In contrast, for the dashed signal, there is destructive interference between components of the signal that reflected from different parts of the rough surface, leading to reduced amplitude and changes in signal shape.

Variations in signal amplitude and shape lead to different thickness measurements. Figure 12(B) shows the Hilbert envelope of the signals in Fig. 12(A) plotted on a dB scale. The 0 dB point on these lines corresponding to the maximum of the Hilbert envelope, which envelope peak detection (EPD) uses to calculate $uLt2$. For the solid signal, this point is at 23.5 mm, compared to 24.5 mm for the dashed signal. Although the actual mean thickness at these two points is about the same, the underlying surface morphology leads to alterations in the reflected signal shape and, therefore, a difference in the measured thickness.

Thickness measurement variation due to changes in pulse shape is determined by the timing algorithm. The black dashed line in Fig. 12 (bottom) is the −10 dB line. The first point where the signal crosses the line, is used by threshold first arrival (TFA) to calculate the thickness. In contrast to EPD, there is only a small difference between the thickness measured from the solid and dashed signals using TFA.

The collected signal can be thought of as superposition of an average coherent pulse, corresponding to the mean thickness of the plate and an incoherent component, introduced by the backscatter from the surface morphology. Averaging signals collected across the plate reveals the shape of the average coherent pulse. As the reflected energy is shared between the incoherent and the coherent components and the fraction of energy of the incoherent component increases with surface roughness, for rougher surfaces one would expect the average pulse to drop in amplitude. This is shown in Fig. 13. For the experimental and simulated data, it can be seen that the average signal amplitude drops rapidly with increasing rms surface variation. The collected signal is being dominated by the random component introduced by the surface roughness. Therefore, for surfaces with larger rms surface variations, there will be much larger changes in signal shape and amplitude, leading to larger variations in the thickness measurement.

One should note that, in Fig. 13, there is some variation in the position of the average signal. This is because the manufacturing process of the rough surfaces is difficult and complex, and offsets between the mean planes of the surfaces could not be avoided. Furthermore, the average pulse shape between the experimental and simulated mean signals (Fig. 13) is different as the transducer used in the experiments was driven by a pulser, as opposed to the Hanning windowed pulse used in the simulations. The bandwidth of both pulses is approximately the same.

### B. ECDFs

Corrosion maps provide an overview of the condition of a component. However, it is difficult to draw quantitative conclusions from them. The CDF, offers a compact presentation of the thickness measurements. Independent statistical distributions can be clearly distinguished, and form the basis for drawing conclusions about the corrosion processes occurring in a vessel.^{5} The CDF can be calculated by sorting the thickness measurements into ascending order and assigning each a rank from 1 to $N$ (where *N* is the number of thickness measurements). The CDF is then

where $X$ is the thickness measurement, $i$ is its rank, and $N$ is the total number of measurements. $F(X)=P(x<X)$ is the probability of measuring a thickness less than $X$. This definition of the CDF has been chosen for consistency with the existing literature on the CDF of thickness measurements.^{5}

Figure 14 shows the ECDFs of the thickness measurements from the three different plates, calculated using EPD. The crosses denote the experimental results and the circles denote the simulated results. There is close agreement between the simulations and the experiments, both showing an increasing spread in the measurements with the increase in rms surface variation. This is expected, as for increasing rms surface variation, the reflected signals will become more incoherent.

The majority of differences between the experimental and simulated ECDFs arises in the tail ($F(x)<10\u22123$), which corresponds to a small number of measurements (0.1%). These measurements correspond to signals which have undergone quite large pulse shape changes, due to the surface roughness and these measurements will be the most susceptible to noise in the experimental set-up.

Figures 15(A) and 15(B), 16(A) and 16(B), and 17(A) and 17(B) show ECDFs for the thickness measurements extracted from the rms = 0.1, 0.2, 0.3 mm plates with different timing algorithms. The graphs labeled A show the experimental results (with crosses) and the graphs labeled B show the simulated results (with circles). Each timing algorithm produces a different tail.

The tail of the distribution can be used for extrapolation purposes.^{5} For example, one could assume that uninspected areas of a component have the same thickness distribution as that measured over a small sample area. The ECDF can be interpreted as the fraction of measurements with less than a given thickness.^{5} From Fig. 17(A), the probability of measuring a thickness less than $x\xaf\u22121\u2009mm$ (at any random measurement point), where $x\xaf$ is the mean thickness, is $F(x=\u22121\u2009mm)\u22480.002$ for EPD and $F(x=\u22121\u2009mm)=0.006$ for cross-correlation (XC). Interpreting these values as a percentage of the area, an inspector would conclude that 0.2% of the structure would have a thickness of less than $x\xaf\u22121\u2009mm$, using EPD, or 0.6% with XC. The actual percentage of the component with less than $x\xaf\u22121\u2009mm$ thickness is 0.05% (from the point cloud).

Clearly, this has consequences for any extrapolation scheme using ultrasonic thickness measurements. First, the minimum thickness in the uninspected area will be underestimated. Second, ultrasonic thickness measurements of the worst case defect for these surfaces will lead to overestimations of the probability of measuring less than a given thickness; the size of the overestimation is determined by the timing algorithm used.

The overestimation will get worse with increasing rms surface variation. In Figs. 15–17 the difference of measured thickness distribution to the point cloud of the actual thickness values, grows with increasing rms surface variation; the rougher the surface, the larger the overestimation of the size of the worst case defect.

### C. Standard deviations of the thickness measurements

ECDFs plotted on a semi-log axis are very good for showing differences in the distribution tails; however, they suppress differences in the bulk of the distribution. For example, it is hard to see quantitatively from Figs. 15–17 how the overall spread in the measurements varies with the choice of timing algorithm or surface roughness. The measured thickness distributions in this paper are all Gaussian, due to the nature of the surfaces. Therefore, calculating the standard deviation of the thickness measurements will give a measure of how the spread in the measurements changes. By examining the standard deviation as a function of surface roughness, one can draw conclusions about how much measurement error could be introduced by each timing algorithm. To determine this, one needs a measure of how much of the standard deviation can be attributed to the surface roughness alone.

The expected standard deviation from the surface roughness can be calculated by considering the patch of surface insonifed by the transducer for each measurement (Sec. III B 2). It is assumed that the thickness which should be measured by the transducer is the mean thickness of this patch. The expected standard deviation of the thickness measurements is the standard deviation of the means of the patches.

In Fig. 18(A) the standard deviation of the thickness measurements is shown as a function of the rms surface variation of each surface, for both the simulated and experimental results. The black dashed line in each figure is the standard deviation expected from the surface roughness. The top graph shows the standard deviation calculated using EPD, the middle using TFA, and the bottom using XC.

For the rms = 0.2 and 0.3 mm plates, the standard deviation of the simulations agrees closely with the experiments, for all of the timing algorithms. However, for EPD and TFA, the standard deviation of the rms = 0.1 mm experimental thickness measurements is slightly higher than the simulated results. This is due to other noise sources in the experimental data which become more significant compared to the signal changes introduced by the roughness when the overall noise is low.

The effect of incoherent noise is less pronounced for XC, as it relies on pulse shape and it is very good at rejecting random noise, while TFA and EPD rely on the signal exceeding the noise floor.^{8} This is shown by the excellent agreement between the simulated and experimental standard deviations for XC in Fig. 18(A). However, it should be noted that thickness measurements extracted using XCs have a larger standard deviation than both TFA and EPD. The increased spread introduced by the use of XC is caused by its reliance on pulse shape. For increasing rms surface variation, the surface roughness has a larger and larger effect on the reflected pulse; this leads to a larger standard deviation.

EPD and TFA perform better than XC with increasing surface roughness, as they are not as susceptible to changes in pulse shape as XC. TFA performs the best out of the three algorithms as the starting point of a signal, is a more stable estimate of the thickness measurement than the peak of a pulse. The peak of the pulse moves around with changes in pulse shape (see Sec. V A). Up to rms = 0.1 mm ($\lambda /12$), the standard deviations of the measurements for all the timing algorithms match up well with the expected standard deviations. However, past this point, the standard deviations increase at a much larger rate than expected; this rate is determined by the choice of timing algorithm. Therefore, one can conclude that up to 0.1 mm ($\lambda /12$) rms surface variation, the spread in the measurements is dominated by the surface roughness, while for rms surface variations greater than 0.1 mm, it is dominated by errors introduced by the timing algorithm.

### D. Frequency dependence

A study of the frequency dependence of the standard deviation of the thickness measurements was also performed, with the same experimental set-up. The plates were scanned using a 3.5 MHz 6 mm diameter longitudinal transducer and the standard deviations of the thickness measurements were calculated. The standard deviations were plotted as a function of the rms surface variation of the plates (Fig. 19). It should be noted that rms surface variation is given here as a fraction of the incident wavelength. The expected standard deviation (black dashed line) was calculated using the field at the backwall for a 5 MHz transducer (Sec. V C); the line for 3.5 MHz has not been included as it does not differ significantly from this line.

It is clear from Fig. 19 that the standard deviation increases as a function of $rms/\lambda $. Up to $rms/\lambda =0.1$ the standard deviation is close to the expected values. However, past this point, it increases faster than anticipated. This suggests that to obtain a thickness distribution that is more representative of the actual surface, one should increase the wavelength of the incident pulse by using a lower frequency transducer. This is at odds with the current paradigm in the NDT industry, where there is a belief that higher frequency always leads to higher accuracy. However, a reduction in frequency may conflict with the need for a higher frequency to resolve the back face echoes in cases where wall-thickness is low.

## VI. CONCLUSIONS

A modified DPSM was used to simulate C-scans of surfaces representative of a corroded engineering component. It was found that the addition of zero pressure boundary points around the outside of the transducer significantly improved the matching of the DPSM solution to an analytical benchmark. Once the DPSM solution was validated by an analytical expression for reflection from flat surfaces, it was used to model the reflections from rough surfaces. Ultrasonic reflections from about 100 000 rough surfaces were modeled and compared to experimental ultrasonic scans of the same surfaces. The statistics of the model and the experiments agreed well, showing that the DPSM simulations are an effective tool to simulate populations of ultrasonic signals from rough surfaces.

The simulated and experimentally acquired ultrasonic signals from surfaces with three different rms surface variations, were analyzed using three different timing algorithms in order to extract thickness measurement data. It was found that the thickness measurement distribution can differ significantly from the actual surface distribution, especially in the tail of the distributions and at larger rms values. Furthermore, the shape of these distributions changed with the choice of the timing algorithm, which implies that the assessment of a component will be dependent on the timing algorithm used to extract the thickness measurements from the A-scans.

The standard deviation of the thickness measurements was also investigated. It was found that up to 0.1 mm ($\lambda /12$) rms surface variation the standard deviation increased in proportion to the change in actual surface roughness. However, for larger rms surface variation the standard deviation increase of measured thicknesses was larger than that of the underlying surface and dependent on the choice of timing algorithm. A study of the frequency dependence of the standard deviation was also performed. It was found that reducing the frequency of the interrogating wave, reduced the error introduced by the timing algorithm and the overall standard deviation of the thickness measurements. This is counterintuitive to the general ultrasonic thinking, where increased frequency is associated with increased resolution and accuracy. The results in this paper show that thickness distributions of rough surfaces might be more precisely assessed when interrogated at lower frequency so that $rms/\lambda <0.1$.

## ACKNOWLEDGMENTS

The authors would like to thank Gordon Davidson (Sonomatic Ltd., Warrington, UK) for his advice and assistance during the course of this project and Simon Burbidge (Imperial College High Performance Computing^{38}) for assistance running the simulations.

## References

*Prevention of Oil Pollution Act*, 1971, http://www.legislation.gov.uk/ukpga/1971/60 (Last viewed 6/27/14).