Two-Dimensional Ultrasonic Flaw Detection ... - Creatis - INSA Lyon

1382

ieee transactions on ultrasonics, ferroelectrics, and frequency control, vol. 44, no. 6, november 1997

Two-Dimensional Ultrasonic Flaw Detection Based on the Wavelet Packet Transform Marc C. Robini, Isabelle E. Magnin, Member, IEEE, Hugues Benoit-Cattin, and Atilla Baskurt Abstract—An important issue in ultrasonic nondestructive evaluation is the detection of flaw echoes in the presence of coherent background noise associated with the microstructure of materials. Many signal processing techniques have proven to be useful for this purpose, but fully 2D flaw detection techniques remain desirable. In this paper, we describe a novel automatic flaw detection method based on the wavelet packet transform, which is particularly well adapted to B-scan image analysis. After a brief review of the essential elements of the theory of wavelets and wavelet packets, a detailed description of the method is provided. The detection process operates on a set of spatially oriented frequency channels, i.e., detail images, obtained from successive wavelet packet decompositions of the initial B-scan. A statistical selection procedure based on the modeling of the detail image histograms retains the useful informationbearing frequency channels. The flaw information is then extracted from these selected channels by means of a specific thresholding scheme. Some experimental detection results in B-scan images of austenitic stainless steel samples comprising artificial flaws are presented.

I. Introduction ltrasonic imaging plays an important role in nondestructive evaluation (NDE) of engineering materials as well as in noninvasive diagnostic medicine. However, the detection capability is often severly limited by the interference noise (i.e., clutter, speckle) produced by the unresolvable scatterers randomly distributed throughout the material. Such coherent interference noise may become significant to the point that it completely masks the target of interest (i.e., flaw, tumor, etc.) in both A-scan signals and B-scan images. As a result, a number of techniques have been proposed to improve the signal-to-noise ratio. Typical approaches include spatial compounding [1]–[3], frequency compounding [4], [5], linear bandpass filtering [6], [7] and split-spectrum processing (SSP) techniques such as averaging, minimization, polarity thresholding and Bayesian detection [8]–[12]. In the NDE field, SSP methods give excellent results but can be sensitive to the choice of parameters [13], [14]. Alternative approaches that show less sensitivity to the environment make use of adaptive filtering, order statistic filtering, and constant false-alarm rate detection [15]–[17]. The use of the wavelet transform

U

Manuscript received December 13, 1995; accepted March 1, 1997. This work was supported by EDF (Electricit´ e de France) under contract #D12160-08 and has been performed in accordance with the scientific trends of the national research group GDR ISIS of the CNRS (Centre National de la Recherche Scientifique). The authors are with CREATIS, INSA 502, 69621 Villeurbanne Cedex, France (e-mail: [email protected]).

was also recently considered [18], [19]. When considering B-scan images, fully 2-D flaw detection techniques are desirable whenever the flaw extends over two or more adjacent A-scans, but it appears that few processing techniques perform simultaneously on both dimensions of the image formed [11], [20]–[22]. In the current work, we propose a novel method for automatic 2-D flaw detection based on the wavelet packet transform. The presence of flaws introduces some modifications of the spectral properties in both directions of the B-scan that can be analyzed separately. In the temporal direction, the physics is governed by the shape of the transmitted sound pulse. If the grain scattering is mainly in the Rayleigh region, the noise spectrum will have most of its power in the high-frequency region of the transducer pass-band [6], [7], [20]. This is not the case for flaw echoes because flaws are generally larger in size than the grain and often behave like geometrical reflectors. In fact, due to the overall filtering effect of attenuation, the flaw spectrum will have a stronger content in the low-frequency region of the transducer pass-band. In the spatial direction, provided that the flaw size is significantly larger than the average grain size of the metal sample being tested, a small shift in the transducer location can result in a significant change of the grain interference pattern, while it has minor effects on the flaw signal [3], [20]. Hence, a relatively flat power spectrum is observed for the time instants corresponding to grain noise only, whereas for the flaw signal, the correlation between the adjacent echoes yields a higher intensity low frequency region [20]. The behavior of the temporal and spatial frequency components of a B-scan image at flaw locations is a productive attribute for detection. Optimal bandpass filtering can provide good results if the information-bearing frequency bands that are dependent on many complex factors (such as the grain size, the flaw size, location and orientation, the scattering properties of the material, the characteristics of the transducer, etc.) are known a priori. Yet, in most cases, this information is unknown and alternative flaw detection techniques such as the ones described in [16] and [17] can be used for A-scan processing. When dealing with a B-scan image, one should clearly consider its horizontal (spatial) and vertical (temporal) directions as preferential. Keeping this in mind, the 2-D separable wavelet transform [23], [24], seems to be an effective tool for B-scan processing. However, the conventional wavelet transform decomposes a signal into a set of frequency channels that have narrower bandwidths in the lower frequency region, which makes it suitable for signals consisting pri-

c 1997 IEEE 0885–3010/97$10.00

robini

et al.:

ultrasonic flaw detection

marily of smooth components. This is not directly applicable to bandpass signals such as ultrasonic RF signals and leads naturally to a generalization of the concept of wavelet bases, namely wavelet packets (WP’s) [25], [26], which provide the required amount of flexibility to analyze B-scan images. The design of the flaw detection algorithm proposed in this paper is motivated by the fact that the description of the B-scan in some particular frequency bands (called the useful information-bearing frequency bands) exhibits a stronger flaw-to-clutter ratio (FCR) than the original signal. As suggested above, this assumption holds if the grain scattering is mainly in the Rayleigh region and if the flaw is larger in size than the grain and is intercepted by several adjacent A-scans. Under these conditions, the received wideband B-scan is first partitioned into several sets of independent, spatially oriented frequency channels (detail images) via successive applications of the WP transform, each of these channels being then a potential candidate for flaw detection. Next, a selection of the useful informationbearing frequency channels, that is the channels with the highest FCR, is performed on the basis of the statistics of the detail images, i.e., their probability density functions (PDFs). Appropriate modeling of these PDFs therefore constitutes a key issue of the method. To complete detection, that is to come up with binary results indicating the location of the flaw echoes that are present in the processing window, a global thresholding scheme is employed to extract the flaw information from the high FCR detail images retained by the selection procedure. However, because of the unimodal characteristic of the detail image PDFs, the common threshold selection methods based on the amplitude features are not adapted. Furthermore, due to the varying statistics of the detail images, the thresholds must be set adaptively. Lastly, a simple combination procedure leads to the final detection results. The paper is organized as follows: in Section II, we briefly review the theory of wavelets and wavelet packets. The application of the WP transform to flaw detection is described in Section III. Section IV presents experimental results for B-scan images of austenitic stainless steel samples (comprising artificial flaws) evaluated with different types of transducers; it appears that the proposed algorithm behaves well in the presence of multiple flaws with different FCR’s in the processing window. Concluding remarks are given in Section V.

1383

[27], [28]:

WTf (a, b) =

+∞ Z f (x)ψa,b (x) dx, −∞

ψa,b (x) = a−1/2 ψ

(1)

x−b . a

In practice, one prefers to write f as a discrete superposition corresponding to a discrete set of continuous wavelets. Of particular interest is the discretization on a dyadic grid (a = 2j , b = 2j k, (j, k) ∈ Z2 ) for which it is possible to construct functions ψ such that the set (ψj,k (x) = 2−j/2 ψ(2−j x−k))(j,k)∈Z2 constitutes an orthonormal basis of square integrable functions over R (i.e., L2 (R)) [27]–[30]. On this discrete grid, the wavelet decomposition of f becomes: XX j dk (f )ψj,k (x) (2) f (x) = j

k

where the wavelet coefficients djk (f ) are inner products of the signal with the wavelet basis functions. djk (f )

+∞ Z = hf, ψj,k i = f (x)ψj,k (x) dx.

(3)

−∞

This discrete wavelet transform has been extensively studied by Mallat [23], [24], who built a complete mathematical theory of multiresolution analysis particularly well adapted to the use of wavelet bases in image analysis. In a multiresolution analysis, one considers the set of vector spaces (Vj )j∈Z . These spaces Vj describe successive approximation spaces · · · V2 ⊂ V1 ⊂ V0 ⊂ V−1 ⊂ V−2 · · · of L2 (R), with resolution 2−j . One also introduces a scaling function φ, together with its dilated and translated versions φj,k (x) = 2−j/2 φ(2−j x − k), (j, k) ∈ Z2 , such that (φj,k (x))k∈Z constitutes an orthonormal basis of the closed subspace Vj . For each j, the ψj,k span a space Wj which is exactly the orthogonal complement in Vj−1 of Vj . Therefore, the wavelet coefficients djk (f ) describe the difference of information between the approximation of f with resolution 2−(j−1) and the coarser approximation with resolution 2−j . To construct the mother wavelet ψ, we may first determine the scaling function φ which satisfies the two scale difference equation: √ X φ(x) = 2 h(n)φ(2x − n) (4) n

II. Wavelet Transform and Wavelet Packets A. A Short Review of Wavelet Analysis The wavelet transform of a function f corresponds to the decomposition of f on the family of wavelets (ψa,b (x))a∈R∗+ ,b∈R generated from one single function ψ (the mother wavelet) by dilatations and translations [23],

where h is the impulse response of a discrete filter defined by h(n) = 1/2hφ(x/2), φ(x − n)i, n ∈ Z. The mother wavelet ψ is related to the scaling function via √ X ψ(x) = 2 g(n)φ(2x − n) (5) n

where g(n) = (−1)n h(1 − n), n ∈ Z. The algorithm proposed by Mallat [23] for the computation of the djk (f ) can

1384


then be summarized by the following two equations: X djk (f ) = g(2k − n)aj−1 n (f ) n

ajk (f ) =

X

h(2k − n)aj−1 n (f )

(6)

n

where the ajk (f ) are coefficients characterizing the projection of f onto Vj . If the function f is given in a sampled form, then one can take these samples for the highest order resolution approximation coefficients a0k and (6) describes a subband decomposition on these sampled values with low-pass filter h and high-pass filter g. As originally investigated by Daubechies [27], it is possible to construct orthonormal wavelet bases in which ψ has finite support, therefore corresponding to FIR filters. Much work investigating the design of such filters can be found in the signal processing literature [31]–[34]. An obvious way to extend the 1-D wavelet transform to higher dimensions is to use separable wavelets [23], [24]. Consider a 1-D scaling function φ(x) and its associated wavelet ψ(x). One can construct four 2-D functions: Φ(x, y) = φ(x)φ(y), Ψ2 (x, y) = ψ(x)φ(y),

Ψ1 (x, y) = φ(x)ψ(y), Ψ3 (x, y) = ψ(x)ψ(y),

(7)

which are orthogonal to each other with respect to integer shifts. The function Φ(x, y) is a separable 2-D scaling function (a low pass filter) and the wavelets Ψi (x, y) are such that the set (2−j Ψi (2−j x − k, 2−j y − 1))i=1,2,3;(j,k,l)∈Z3 is an orthonormal basis of L2 (R2 ). The approximation of a signal f (x, y) at resolution 2−j is characterized by the set of inner products Aj f = hf (x, y), 2−j Φ(2−j x − k, 2−j y − l)i (k,l)∈Z2 , (8) and the difference of information between Aj−1 f and the coarser approximation Aj f is given by the detail images Dij f = hf (x, y), 2−j Ψi (2−j x − k, 2−j y − l)i (k,l)∈Z2 , i = 1, 2, 3. (9) This solution corresponds to a separable 2-D filter bank with subsampling by two in each dimension. Fig. 1 represents one stage in such a decomposition of an image, together with the corresponding partition of the frequency plane. Starting with the function f given in a sampled form (A0 f ), this decomposition scheme is recursively performed to the output giving the low resolution subimage Aj f . It leads to the conventional wavelet transform which can be interpreted as an octave band signal decomposition. B. Wavelet Packets WP’s were introduced by Coifman et al. [25] and Wickerhauser [26], as a family of orthonormal bases for discrete functions of RN , including the wavelet basis and the shorttime-Fourier-transform (STFT)-like basis as its members.

Fig. 1. (a) Separable 2D filter bank corresponding to a separable wavelet basis, the filters h and g are, respectively, a half-band lowpass filter and a half-band highpass filter. Aj−1 f is the initial image corresponding to resolution 2−(j−1) , Aj f is the low resolution subimage (2−j ), and (Dij )i=1,2,3 are the detail images corresponding to the information visible at resolution 2−(j−1) . The corresponding partition of the frequency plane is indicated in (b).

WP’s represent a generalization of the method of multiresolution decomposition in the sense that the WP basis functions (ϑm )m∈N can be generated from a given function ϑ0 according to: √ X h(n)ϑm (2x − n) ϑ2m (x) = 2 n

ϑ2m+1 (x) =

√ X 2 g(n)ϑm (2x − n)

(10)

n

where the function ϑ0 can be identified with the scaling function φ and ϑ1 with the mother wavelet ψ. Then, the library of WP bases can be defined to be the collection of orthonormal bases of L2 (R2 ) composed of functions of the form ϑj,k,m (x) = 2−j/2 ϑm (2−j x − k), (j, k, m) ∈ Z2 × N. Each element of the library is determined by a subset of the indexes j, k, and m (for instance, a conventional wavelet basis corresponds to the collection of indexes (j, k, 1), (j, k) ∈ Z2 ). The set (ϑj,k,m (x))k∈Z is an orthonormal basis of a subspace Vj,m of L2 (R) such that Vj,2m+1 is the orthogonal complement of Vj,2m in Vj−1,m , where Vj,0

robini

et al.:


1385

can be identified with the closed subspace Vj and Vj,1 with Wj . By analogy with (7), one can construct the set of functions (Vm,n (x, y) = ϑm (x)ϑn (y))(m,n)∈N2 ,

(11)

where V0,0 = Φ, V0,1 = Ψ1 , V1,0 = Ψ2 , and V1,1 = Ψ3 stand as particular cases. A 2-D WP basis is then composed of functions of the form:

TABLE I 8-taps Daubechies Filter Coefficients. h(0) h(1) h(2) h(3) h(4) h(5) h(6) h(7)

−0.0535744507090 −0.0209554825625 0.3518695343280 0.5683291217040 0.2106172671020 −0.0701588120895 −0.0089123507210 0.0227851729480

Vj,k,l,m,n (x, y) = 2−j Vm,n (2−j x − k, 2−j y − 1) (12) where the collection of indexes (j, k, l, m, n) ∈ Z3 × −j 2 that the intervals 2 m, 2−j (m + 1) × N −j is such −j 2 n, 2 (n + 1) form a disjoint cover of [0, +∞) × [0, +∞) and k, l range over all the integers. Hence, any triplet (j, m, n), (m, n) 6= (0, 0) thus defined gives rise to a detail image: j Dm,n f = (hf (x, y), Vj,k,l,m,n (x, y)i)(k,l)∈Z 2

(13) which describes f (x, y) in the frequency bands [−(m + 1)2−j πx , −m2−j πx ] × [−(n + 1)2−j πy , −n2−j πy ] ∪ [m2−j πx , (m + 1)2−j πx ] × [n2−j πy , (n + 1)2−j πy ]. From a practical point of view, the key difference between the conventional wavelet transform and the WP transform is that the decomposition scheme [Fig. 1(a)] is no longer simply applied to the output giving the low frequency subimage. Instead, it can be applied to any output, thus leading to a quadtree structure representation. Let us notice finally that STFT-like bases correspond to the particular case of a regular tree, i.e., when all subimages are decomposed in each scale to achieve full decomposition, or equivalently, uniform frequency resolution.

III. Wavelet Packets Based Ultrasonic Flaw Detection Due to the separability of the B-scan into orthogonal directions and based on the obervation that its partial spectra (the temporal and spatial frequency components) behave in a distinct manner at flaw locations [3], [6], [7], [20], the 2-D separable wavelet transform appears to be wellsuited to B-scan image analysis. The conventional wavelet transform [23], [24], recursively decomposes subsignals in the low frequency channels. However, because the information contained by ultrasonic RF signals is located in the middle frequency channels, further decomposition just in the lower frequency region does not help much for the purpose of flaw detection. Thus, an appropriate way to analyze B-scan images is to allow the decomposition of any frequency channel such as the WP transform [25], [26]

does. As we are concerned with detection, redundancy is not prejudicial, it is even desirable. Therefore, we choose to extract the flaw information from the set of detail imj f )j=1,...J;m,n=0,...2j −1,(m,n)6=(0,0) resulting from ages (Dm,n the decompositions of the initial B-scan A0 f into STFTlike bases of successive depths 1, . . . J (i.e., resolutions 2−1 , . . . 2−J ). From a practical point of view, the size of the smallest subimages at resolution 2−J may be used as a stopping criterion for further decomposition. But in order to limit the performance decrease introduced by the resolution loss, we require the numerical period tc (expressed in pixels) associated with the nominal center frequency fc of the transducer to be at least one pixel wide at the smallest resolution. Therefore, if fs denotes the sampling rate for the A-scan signals that constitute the B-scan, then tc = fs /fc has to be such that 2−J tc ≥ 1 and the stopping criterion is given by: log2 (fs /fc ) − 1 < J ≤ log2 (fs /fc ).

(14)

An example of wavelet packet decompositions of a Bscan into STFT-like bases of depths 1 and 2 is shown in Fig. 2. We used the 8-taps Daubechies filter associated with the “least asymmetric” compactly supported wavelet [27] which normalized coefficients are listed in Table I. The data were collected from an austenitic stainless steel sample with 2 mm diameter holes representing flaws, using a broad-band, focused angle probe (see Section IV for more details: Table III, B-scan #5). As mentioned earlier, the proposed flaw detection method operates in four steps. A simplified diagram of the process is provided in Fig. 3. The selection of the useful information-bearing frequency channels, to be described in Section III.B, is based on the modeling of the detail image PDFs by generalized Gaussian functions, which is discussed in Section III.A. The adaptive thresholding scheme and the final combination procedure are laid-out in Sections III.C and III.D, respectively. A. Statistical Properties of Wavelet Coefficients Because the pixels of the detail images are the decomposition coefficients of the initial image in an orthonormal family, their distribution could exhibit any shape. However, it can be verified in practice that the PDF of the detail images resulting from the decomposition of a B-scan

1386


Fig. 2. WP decompositions of a B-scan image (a) (see Section IV and Table III (B-scan #5) for experimental description) into STFT-like bases of depths 1 (b) and 2 (c) (resolutions 2−1 and 2−2 ). According to the scheme presented in Fig. 1(a), the B-scan image is first decomposed into A1 f , D11 f , D21 f , and D31 f , respectively, located in the upper left, upper right, lower left, and lower right quadrant of (b). Each of these four subimages is then decomposed in the same way to obtain the full decomposition of depth 2 and so on. The wavelet coefficients in the detail images are displayed in absolute value, the useful-information-bearing frequency bands appear clearly.

Fig. 4. (a) High FCR detail image resulting from the WP decomposition of depth 2 of the B-scan in Fig. 2. (b) Associated normalized histogram p(x) and its approximation model p˜(x) (20) (α = 0.047, β = 0.830).

Fig. 3. Principle of the flaw detection method.

in a WP basis are symmetrical peaks centered in zero. Indeed, these PDF can be modeled by generalized Gaussian functions, as was already found experimentally in [23] and [35] for other types of signals. The generalized Gaussian law is given by:

p˜(x) =

β β e−(|x|/α) 2αΓ(1/β)

(15)

where Γ(.) is the standard Gamma function. This is a twosided symmetric density with two distributional parameters α and β, such that α controls the variance and β modifies the decreasing rate of the central peak. The general formula (15) contains two interesting examples as particular cases: β = 2 leads to the well-known Gaussian PDF and β = 1 leads to a Laplacian PDF. The parameters α and β are computed from the mean absolute deviation m1 and the variance m2 of the detail image: +∞ Z m1 = |x|p(x) dx, −∞

+∞ Z m2 = x2 p(x) dx. −∞

(16)

robini

et al.:


1387 where sign(.) is the signum function and P (.) denotes the incomplete Gamma function,

P (a, x) =

1 Γ(a)

Zx

e−ν ν a−1 dν,

a > 0.

(23)

0

Fig. 5. (a) Grain only detail image resulting from the WP decomposition of depth 2 of the B-scan in Fig. 2. (b) Associated normalized histogram p(x) and its approximation model p˜(x) (20) (α = 0.289, β = 1.842).

By replacing the PDF p(x) of the detail image by (15) and changing variables in these two integrals according to: β α 1−β x ⇔ dx = u β du, α β

u=

(17)

one can derive that m1 = α

Γ(2/β) Γ(1/β)

and m2 = α2

Γ(3/β) . Γ(1/β)

(18)

Thus, β = G−1

2 m1 , m2

where G(x) =

Γ(2/x)2 Γ(3/x)Γ(1/x)

(19)

and α follows directly from either of the two equations in (18). Typical examples of PDF modeling are shown in Figs. 4 and 5. The quality of the approximation provided by the generalized Gaussian model can be assessed with the help of the measure: V =

max

(F (x) − F˜ (x)) +

Decompositions in STFT-like bases of depths 1, 2, and 3 (that is resolutions 2−1 , 2−2 , and 2−3 , respectively) were applied on six flaw-plus-grain B-scan images obtained with different transducers, using a 8-taps Daubechies wavelet filter (Table I). The PDFs of the resulting detail images were modeled by generalized Gaussian functions (15), (18), (19) and the corresponding distance measures (20) were then computed and averaged for each resolution level. We respectively obtained E(V ) = 0.018, 0.030 and 0.039 for resolutions 2−1 , 2−2 , and 2−3 , that is 1.8%, 3%, and 3.9% of the dynamic range of the CDFs. These results clearly demonstrate that the generalized Gaussian PDF model (15) is well adapted to the detail image distributions. Also, as suggested by the two examples considered in Figs. 4 and 5, a faster decreasing rate is observed for the PDFs of the detail images with high FCR, which translates to smaller values of the estimated β parameters. To demonstrate this, the high FCR detail images were manually selected and the associated mean β value was found to be 0.770, whereas E(β) = 1.644 for the remaining low FCR or grain only subimages; meaning, therefore, that the use of β can be very helpful for the selection of the useful information-bearing frequency bands. B. Selection of the Detail Images The contents of a detail image with high FCR primarily consists of a majority of low amplitude wavelet coefficients associated with clutter and a few coefficients with significantly higher amplitudes representing flaws. Hence, as shown above, the corresponding PDF exhibits an important decreasing rate such that, for some fixed threshold Tβ > 0, β ≤ Tβ can be used as a first selection criterion. Let us also introduce the measure:

−∞