An H∞ deconvolution filter and its application to ultrasonic nondestructive evaluation of materials

ISA Transactions 38 (1999) 323±335

www.elsevier.com/locate/isatrans

An H1 deconvolution ®lter and its application to ultrasonic nondestructive evaluation of materials Timothy C. Hanshaw a, Michael J. Anderson b, Chin S. Hsu a,* a

School of Electrical Engineering and Computer Science, Washington State University, Pullman, WA 99164-2752, USA b Department of Mechanical Engineering, University of Idaho, Moscow, ID 83844-0902, USA

Abstract Deconvolution can be a valuable technique for interpreting results of ultrasonic nondestructive evaluation (NDE) tests of materials. This is especially true for state of the art hybrid materials. In this paper, a new H1 deconvolution ®lter is presented. The ®lter is applied to typical ultrasonic NDE data, including through-transmission data for aluminum and composite samples. The results are compared to those obtained from a Wiener deconvolution ®lter. The performance of the H1 ®lter is as good or better than the performance of the Wiener ®lter. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Bounded real lemma; H1 deconvolution ®lter; Nondestructive evaluation of materials

1. Introduction It is often necessary to evaluate the properties of a material sample without a€ecting the sample's usefulness. This process is called Non-Destructive Evaluation (NDE) of materials. NDE can be particularly important in the evaluation of layered metal and composite materials. Since the manufacturing processes used to make these materials are still under development, the ability to evaluate the quality of the part prior to its use and over the course of its life is critical. One common method of performing a nondestructive evaluation is to apply an ultrasonic acoustic pressure pulse to the sample in question; the resulting pressure pulses re¯ected from or * Corresponding author. Tel.: +1-509-335-2342; fax: +1509-335-3818. E-mail addresses: [email protected] (C.S. Hsu), thanshaw @eecs.wsu.edu (T.C. Hanshaw), [email protected] (M.J. Anderson).

transmitted through the sample are then evaluated to determine the properties of the sample [1]. Ultrasonic NDE techniques do, however, su€er a generic limitation in resolution. In general, it becomes dicult to identify a defect if the defect size is smaller than the width or wavelength of the probing ultrasonic beam. This limitation extends to the testing of thin layered materials. Echoes will re¯ect from both sides of the thin layer; these echoes tend to overlap unless the wavelength of the beam is smaller than the thickness of the layer. The individual layers in the composite materials of interest in this paper can be very thin. As an example, 5 mil thick titanium layers are being considered for use in these materials; in order to have the wavelength smaller than the layer thickness, a frequency of 50 MHz would be required. Unfortunately, high frequency pulses attenuate rapidly when propagating through a sample; it is therefore dicult to identify defects deep within a sample when using these high frequencies. Obviously, it is necessary to use more sophisticated signal processing

0019-0578/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S0019-0578(99)00034-8

324

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

techniques to extract the desired information from these samples, while using lower-frequency inputs. Deconvolution methods have received considerable attention as one approach to the problem of extracting information about defects from ultrasonic NDE test data, particularly when the defect size is of the same order as the wavelength of the ultrasonic beam. In essence, deconvolution is the process of obtaining an estimate of the input to a system, given the measured output from the system and a model of the system. In the current application, this ``input'' will consist of the impulse response of the test sample and the ``system'' is the transducers and associated electronics being used to apply and measure the ultrasonic pulse. Refs. [2] and [3] provide overviews and comparisons of several deconvolution methods and their application to NDE problems. Ref. [4] discusses the application of deconvolution to another layered system, namely the earth's crust. These methods typically involve the minimization of an L1 or L2 norm of the error in the input estimate. The L2 norm is essentially the total energy in the signal. The L1 norm is computed by integrating the absolute value of the function. H1 deconvolution methods have never, to the authors' knowledge, been applied to these types of problems. The H1 norm based methods require that the ratio of the L2 norm of the measured output to the L2 norm of the system inputs be bounded by some scalar quantity. H1 type approaches have several advantages over L2 type methods: the H1 methods require no limitations or assumptions about the properties of the input or the noise, the input pressure pulse is not required to be minimum phase, and robust implementations of H1 type ®lters are available. The H1 type ®lters are not as computationally expensive as the L1 norm methods. One distinct disadvantage of the H1 approaches is that they do not perform optimally in the presence of white gaussian measurement noise. The goal of this paper is to apply H1 deconvolution methods to the problem of nondestructive evaluation of layered composite materials in order to assess their performance in this important class of problems. Aluminum samples of various thicknesses were tested at several pulse frequencies, in

order to determine the resolution achieved by the H1 methods. A composite material, consisting of a 16.2 ml thick graphite ®ber reinforced epoxy (GFRE) layer sandwiched between two 5 ml thick titanium layers, was also tested. The overall organization of this paper is as follows. A brief outline of the general testing approach and an overview of the basic analysis techniques to be used are provided in Section 2. This is followed by a discussion of the H1 deconvolution ®lter to be used in this paper (Section 3). Results obtained with the H1 ®lter are presented in Section 4 and compared with results obtained with a Wiener deconvolution ®lter, for the aluminum and composite (titanium±graphite±titanium) samples. Conclusions based on the current analysis and recommendations for future work are provided in Section 5. A detailed derivation of the H1 deconvolution ®lter used here is provided in the Appendix. 2. Problem formulation This section provides a brief overview of the methods used to acquire the test data and de®nition of basic variables. The overall approaches toward signal processing and deconvolution are presented in the context of the test situation. Finally, the method used to separate the transducer characteristics from the characteristics of the sample being tested is presented. The test procedure analyzed here consists of immersing the sample to be examined in a water bath. An ultrasonic acoustic pressure pulse is applied to one side of the sample. This pressure input results in pressure waves being transmitted through and re¯ected from the sample. These transmitted and re¯ected pressures are measured and analyzed to determine the sample's properties. The re¯ected pressure pulse is commonly referred to as the ``pulse-echo'' pressure and the transmitted pressure is referred to as the ``throughtransmission'' pressure. A block diagram of the overall test set up, displaying the primitive variables of interest from a signal processing standpoint, is depicted in Fig. 1. For signal processing purposes, the experimental set up of Fig. 1 can be thought of as a series of impulse responses of the

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

various system components. Convolving these impulse responses together results in a mathematical description of the measured output of the system of Fig. 1, as follows: y…t† ˆ …t† h…t† gi …t† go …t† ‡ w…t†

…1†

where y(t) is the measured system output (the recorded voltages), h(t) is the sample's impulse response, gi(t) is the pressure input system impulse response, go(t) is the measurement system impulse response, and w(t) is the measurement noise. It can be seen from Eq. (1) that the measured response voltages are assumed to be the result of a unit impulse, (t). This is for mathematical convenience; no impulse is physically generated in the system. The deconvolution problem can now be formulated as shown in Fig. 2. The object of the deconvolution ®lter is to remove the e€ects of the measurement noise [w(t)] and the known system [gi(t)*go(t)] from the measured signal y(t) in order to obtain an estimate of the system input, uÃ(t).The goal is to minimize the estimate error, the di€erence

325

between u(t) and uÃ(t), in some sense. Since the assumed input to the system is an impulse, estimating the system input will be equivalent to estimating the impulse response of the sample, h(t). This is the desired re¯ectivity or transmittivity sequence of the sample. The deconvolution approach outlined in Fig. 2 requires an impulse response which does not include the sample being tested (the [gi(t)*go(t) impulse response block in Fig. 2). In this paper, this impulse response is obtained via testing. A ``reference'' specimen, with impulse response x(t), is tested. This specimen can be no sample at all, a sample with known properties, or a sample deemed to be ``¯aw free''. The data acquired with this reference specimen can be written as: yref …t† ˆ …t† gi …t† x…t† go …t† ‡ wref …t†

…2†

System identi®cation techniques can be used to obtain the impulse response gi(t)*x(t)*go(t) from the measured output.

Fig. 1. Experimental test set-up block diagram.

Fig. 2. Deconvolution block diagram.

326

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

Now let the impulse response of our test sample of interest, h(t), be written as: …3† h…t† ˆ x…t†
…t† The impulse response
(t) is relative to the reference specimen.
(t) may or may not have any physical meaning, depending on the choice of the reference specimen. Substituting Eq. (3) into (1), the test sample output can be written as: y…t† ˆ …t† gi …t† x…t†
…t† go …t† ‡ w…t† ˆ href …t†
…t† ‡ w…t†

…4†

where href(t)=gi(t)*x(t)*go(t) is the impulse response of the instrumentation system, including the reference specimen. The system has now been cast in a form in which the impulse response href(t) is known from testing on a reference specimen. The deconvolution block diagram for the system being analyzed becomes, therefore, as shown in Fig. 3. A summary of the basic approach used here toward applying deconvolution methods to determine the re¯ectivity or transmittivity sequence of a material can now be outlined as follows: 1. Acquire ultrasonic test data on a reference specimen, yref(t). 2. Use the reference specimen data from above to generate a model of the system impulse response [href(t), in Fig. 3], using system identi®cation techniques. 3. Acquire ultrasonic test data on the test sample whose properties are to be investigated [y(t) of Eq. (1)]. 4. The deconvolution ®lter is applied as shown in Fig. 3, in order to obtain an estimate of the input to the test sample, uÃ(t), given the sample's output, y(t), and the reference wavelet impulse

response, href(t). This input estimate corresponds to an estimate of the impulse response of the test sample,
(t), with respect to the reference sample and can be interpreted accordingly. 3. H1 deconvolution ®lter The discussion of of the previous section was for a physical, continuous-time system. The data acquisition system will, in general, sample the continuous-time data. This will result in only discrete-time data being available for analysis purposes. For that reason, all analysis methods presented here are for discrete-time data. This section provides a brief overview of the H1 deconvolution ®lter used in the subsequent examples. A detailed derivation of the H1 ®lter parameters is provided in the Appendix. Consider a discrete, time-varying system, described by the state space model, x…k ‡ 1† ˆ A…k†x…k† ‡ B1 …k†u…k† ‡ B2 …k†w…k†; x…0† ˆ x0

…5†

y…k† ˆ C…k†x…k† ‡ D…k†w…k† where k2[0,Nÿ1] is an integer index, x2Rn is the state, u2Rm is the input, y2Rp is the measured output, w2Rq is the disturbance input, and x02Rn is the initial state of the system. The matrices A(k), B1(k), B2(k), C(k) and D(k) are bounded functions of k and have appropriate dimensions. Note that the output y is assumed to have no direct information about the input to the system. The disturbance input, w(k), is assumed to have bounded energy in L2[0,Nÿ1], i.e. ÿ
1=2 k w k2 ˆ wT …k†w…k† < 1: The form of the deconvolution ®lter chosen here is

Fig. 3. Deconvolution block diagram.

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

x^ …k ‡ 1† ˆ A^ …k†x^ …k† ‡ B^ 1 …k†u0 …k† ‡ B^ 2 …k†y…k ‡ 1†; x^ …0† ˆ x^ 0 u^ …k† ˆ C^ …k†x^ …k† ‡ D^ 1 …k†uo …k† ‡ D^ 2 …k†y…k ‡ 1† …6† n

where xÃ2R is the estimated state of the system, uÃ2Rm is the estimate of the input to the system, and AÃ(k), BÃ1(k), CÃ(k), DÃ(k) and DÃ2(k) are the deconvolution ®lter parameters to be determined. u0(k) is any known portion of the input; this assumes that the input can be written in the form u…k† ˆ uo …k† ‡ uu …k†

…7†

uu(k) is the unknown portion of the input, and is of primary interest in this paper. According to the ®lter structure of Eqs. (6), the estimate of the system's state is dependent upon the state estimate at the previous time step and the output at the current time step. The estimate of the input to the system is dependent upon the current state estimate and the output at the next time step. In order to estimate the input at a given time step, then, one requires information about the current state of the system, as well as information as to where the system is being driven. This would be the expected form, given Eqs. (5); the actual state of the system at the next time step is dependent upon the previous state and the previous input. The deconvolution ®lter structure of Eqs. (6) is therefore reasonable from a physical standpoint. One drawback to the ®lter structure is its lack of ``smoothing'' Ð only the current and previous time steps are used in the state estimates. It is obvious that the presence of white noise in the output will severely disrupt the ®lter's performance, especially since the ®lter is provided with no information about the noise statistics. The use of smoothing techniques [5,6] have the potential for signi®cantly improving the ®lter's performance. The performance measure for the ®lter of Eqs. (6) is k u^ ÿ u k22 < 2 SUP w;uu 2L2 ‰0;Nÿ1Š k w k22 ‡ k uu k22

…8†

for x0 known. is any pre-speci®ed positive number. This performance measure requires that the 2-norm/2-norm system gain Ð which is, in

327

essence, the H1 norm of the system [7] Ð be less than a pre-speci®ed positive value, for any choice of input to the system. The system is de®ned from the disturbance and unknown inputs to the estimation error. The only requirement is that the input have bounded energy. Since we always deal with ®nite-duration signals, this requirement is automatically met. The solution to the H1 deconvolution ®ltering problem can be obtained by applying the bounded real lemma for time-varying systems [8]. Application of this approach to the ®lter de®ned in Eq. (6) results in the ®lter parameters shown in Eq. (9). A derivation of these parameters is provided in the Appendix. The ®lter shown here allows for timevarying systems; the current work does not take advantage of this capability, although time-varying and spatially-varying systems do appear in NDE applications [2]. A^ …k† C^ …k† D^ 1 …k† B^ 1 …k† B^ 2 …k† D^ 2 …k†

A…k† ÿ B^ 2 …k†C…k ‡ 1†A…k† ÿD^ 2 …k†C…k ‡ 1†A…k† Im ÿ D^ 2 …k†C…k ‡ 1†B1 …k† ^ ÿB1 …k† ÿ B2 …kT†C…k T‡ 1†B1 …k† A…k†Q…k†A …
k†C …k ‡ 1† ‡BA …k†DTA …k† Rÿ1 …k† ˆ BT1 …k†CT …k†Rÿ1 …k†

ˆ ˆ ˆ ˆ ˆ

…9†

where BA …k† DA …k† R…k†

ˆ ˆ ˆ

Q…k ‡ 1†

ˆ

B…k†

ˆ

‰B…k† 0Š ‰C…k ‡ 1†B…k†D…k ‡ 1†Š DA …k†DTA …k† ‡ C…k ‡ 1† A…k†Q…k†AT …k†CT …k ‡ 1† A…k†Q…k†AT …ÿk† ‡ BA …k†
T T T Tÿ1 1 BA …k† ÿ A…k†Q…k†A …k†C …k ‡ 1 T ÿ1 ‡BA …k†Tÿ1 † 1 DA …k††T2 …k†…C…k ‡ 1
T T A…k†Q…k†A …k† ‡ DA …k†Tÿ1 1 BA …k† ‰B1 …k† B2 …k†Š …10†

with Q(0)=0 for the case of known initial conditions (i.e. x^ 0 ˆ x0 ). T1 and T2(k) are invertible matrices, where T1 ˆ Im‡2q ÿ ÿ2 LT L T T T2 …k† ˆ DA …k†Tÿ1 1 DA …k† ‡ C…k ‡ 1†A…k†Q…k†A …k†

CT …k ‡ 1† L ˆ ‰Im 0 0Š …11†

328

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

Fig. 4. H1 deconvolution ®lter structure.

A block diagram illustrating the ®lter structure is shown in Fig. 4.

Table 1 Sample properties and test conditions discussed in this paper Material

4. Results and discussion This section provides deconvolution results using the H1 ®lter discussed previously. The H1 ®lter results are compared with results obtained with a Wiener deconvolution ®lter. Stringent metrics to quantitatively assess the performance of the ®lters presented here are currently being formulated. Assessments of the ®lters' performance provided in this section are primarily subjective, although comparisons of ®lter results with independently calculated physical parameters are made where feasible. For the test data, the state space model of the input pressure pulse required by the H1 ®lter must be determined from the available data, per the previous discussion. In this paper, the reference signal is always recorded without a sample present. Identi®cation of the appropriate statespace model was performed using the eigensystem realization algorithm (ERA) [10±12]. To date aluminum, titanium, and composite (titanium±graphite ®ber reinforced epoxy) materials have been tested. The composite material was supplied by the Boeing Aircraft Company; material properties shown in Table 1 were provided by the Boeing company. Samples analyzed in this

Thickness Sound (ml) speed (m/s)

Aluminum 77.5 Aluminum 16 Composite 5±Ti 16.2±Gr

6240 6240 6100±Ti 2580±Gr

Density (kg/m3)

Input `/l center frequency (MHz)

2700 25 2700 15 4711±Ti 25 1589±Gr

7.67 0.95 0.58±Ti 4.28±Gr

paper include 77.5 and 16 mil aluminum and a composite made of a 16.2 ml GFRE (graphite ®ber reinforced epoxy) layer sandwiched between two 5 ml titanium layers. Data were taken with single-cycle tone bursts of frequencies 15 and 25 MHz supplied to the transducers. The frequency content of the waveform is determined primarily by the transducer characteristics. Throughout this section, the referenced ``nominal input frequency'' refers to the tone burst frequency. Table 1 shows the test samples and conditions which were analyzed for this paper, along with the ratio of the sample's thickness to the wavelength of the input pulse center frequency; this parameter provides an idea of the amount of resolution a pulse can provide over a layer. Each of the test cases outlined in Table 1 is discussed in the following subsections.

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

4.1. 77.5 ml Aluminum, 25 MHz input center frequency The ®rst case discussed is for 77.5 ml thick aluminum. The nominal center frequency of the input is 25 MHz. For this case, the `=l of the sample is 7.67. Re¯ections from within the sample will therefore be spaced far apart in time relative to the duration of the pulse and will not overlap. The resulting through transmission pressure will be easily interpretable from an intuitive standpoint. This case will provide a valuable benchmark for testing the system identi®cation and deconvolution algorithms to be used, before their application to cases where internally re¯ected waves overlap and obscure the physical meaning of the data. Fig. 5 shows the reference (no sample) voltage wave form and the corresponding waveform with the sample in place. The sampling time step size is 1.25 ns; the Nyquist cuto€ frequency is 400 MHz. For the H1 ®ltering case, the ®rst step is to generate a state space model of the reference wavelet. Discrete Fourier Transforms of the reference data indicate that there is essentially no input at frequencies above about 30 MHz, therefore both the reference and through-sample data were low-pass ®ltered using a 10th order Butterworth ®lter; the cuto€ frequency of the ®lter was selected

Fig. 5. Reference signal (top) and through transmitted signal (bottom) for 77.5 ml aluminum sample, 25 MHz (nominal) input center frequency.

329

to be 30 MHz (7.5% of the Nyquist frequency). The ERA was then used to determine a state space model of the reference (no-sample) wavelet. A 10th order model was used to approximate the reference signal; Fig. 6 shows the impulse response of the model vs. the ®ltered reference data. The agreement between the two is excellent and is typical of the level of agreement for all studies conducted. It should be noted that the model is non-minimum phase (the model has zeros outside the unit circle); this is also typical of the models generated in this research. The ®ltered output data were then processed using the H1 deconvolution ®lter. Two representative transmission sequences thus obtained are shown in Fig. 7. Increasing the ``D'' term results in improved signal-to-noise characteristics in the deconvolved signal, at the expense of reduced overall amplitudes. For the Wiener ®lter analysis, the only preprocessing of the data was to remove some of the leading ``zeros'' from the reference signal. This was done in order to avoid shifting the initial re¯ection impulse to a negative time, since the wave propagation speed through the sample is higher than that through water. This e€ect is readily apparent from Fig. 5; the ®rst incidence of a pulse comes earlier in the ``sample'' data than in the ``no-sample'' data. In order to keep the time series lengths the same for the two cases, a similar amount of data was discarded from the end of the ``sample'' data. Fig. 8 shows the deconvolved

Fig. 6. Impulse response of identi®ed 10th-order state space model vs no-sample (reference) data. (Data = +). 77.5 mil Al, 25 MHz case.

330

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

Fig. 7. Deconvolved transmission sequence from H1 ®lter. 77.5 ml Aluminum sample, 25 MHz (nominal) input center frequency: (a) D=0.5, =10; (b) D=2.0, =10.

4.2. 16 ml Aluminum, 15 MHz input center frequency

Fig. 8. Deconvolved transmission sequence from Wiener ®lter. 77.5 ml Aluminum sample, 25 MHz (nominal) input center frequency.

transmission sequence as determined by the Wiener ®lter. The ``sinc'' type behavior of each pulse is readily apparent. Qualitatively, the deconvolved sequences of Figs. 7 and 8 show some improvement over the raw signal (Fig. 5), in the sense that the initial portion of the pulse is more strongly emphasized in the deconvolved responses. Overall, the deconvolved pulses exhibit a more impulse-like behavior than the raw data, as desired.

The next case consists of a 16 ml thick aluminum sample. The (nominal) input pulse center frequency is 15 MHz; this combination results in an `=l of 0.95. For this case, internal re¯ections from the sample overlap, obscuring the meaning of the output pressure measurement. Fig. 9 shows the reference wave form (taken with no sample) and the corresponding throughtransmission wave form for the 16 mil aluminum sample. The raw data obviously does not provide much direct information about the sample characteristics. For this case, the sampling time step size is 0.5 ns; the Nyquist frequency is 1000 MHz. As in the previous case, the ERA was used to generate a state space model of the reference wavelet. Prior to using the ERA, the data were low-pass ®ltered with a 12th order Butterworth ®lter with cuto€ frequency of 35 MHz and decimated by a factor of 2. Without decimation, the size of the Hankel matrices used by the ERA became excessively large and performance of the method was poor. A 10th order model was used to replicate the reference pulse. Fig. 10(a) shows the deconvolved transmission sequences obtained from the H1 ®lter, for the ®rst 1000 data points. Peaks are readily apparent on Fig. 10(a), indicating

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

the location of strong transmission coecients. These peaks are approximately 120 time steps (120 ns) apart. The expected (2-way) wave travel time through the sample is 126 ns, based on an acoustic velocity in the sample of 6420 m/s. For the Wiener ®lter analysis, the only preprocessing of the data was to remove leading ``zeros'' from the input data, as was done on the previous case. The through-transmission data was truncated to keep the time series lengths the same for both cases. Fig. 10(b) shows the deconvolved transmission

331

sequence obtained from the Wiener ®lter. The peaks on the Wiener ®lter results are approximately 240 time steps (120 ns) apart. This agrees well with both the H1 ®lter results and the expected value. Results obtained with the H1 ®lter and the Wiener ®lter are comparable, as can be seen from Fig. 10. The primary di€erence is that the H1 ®lter results exhibit two distinct oscillations between primary peaks, while any oscillations between peaks in the Wiener ®lter results appear to have coalesced. It is unclear whether one approach provides a distinct advantage over the other. 4.3. Ti-GFRE-Ti composite, 25 MHz center frequency

Fig. 9. Reference signal (top) and through-transmitted signal (bottom) for 16 ml aluminum sample, 15 MHz (nominal) input center frequency.

The ®nal case examined consisted of a 16.2 mil thick GFRE layer sandwiched between two 5 ml thick titanium layers. The center frequency of the input pulse was 25 MHz (nominal). Fig. 11 shows the reference (no sample) wave form and the corresponding waveform with the sample in place. The sampling time step size is 1.25 ns and the Nyquist frequency is 400 MHz. The procedure for performing the H1 analysis is as usual. The data was decimated by a factor of two and a 12th order Butterworth ®lter, with cuto€ frequency of 28 MHz, was applied to the data. A

Fig. 10. Deconvolved transmission sequences for 16 mil aluminum, 15 MHz input frequency, H1 ®lter data decimated2: (a) H1 ®lter, D=0.5; =10. (b) Wiener ®lter.

332

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

10th order model of the ®ltered reference waveform was generated using the ERA. The deconvolved transmission sequences obtained from the H1 and Wiener ®lters are shown in Fig. 12. It is dicult to interpret the data of Fig. 12, as its behavior is rather complex. The primary signature of the transmission sequence is an oscillation with a period of approximately 42 ns which are presumably due to the wave re¯ections in the titanium layers; the 2-way wave travel time in the 5 ml titanium layers is estimated as 41.6 ns. Secondary oscillations are also apparent, with a period of

approximately 120 time steps (300 ns). This is most likely the response of the ®lter to the graphite layer Ð the calculated 2-way travel time in the graphite is approximately 319 ns. None of the above observations is signi®cantly more obvious from the deconvolved transmission sequences than they are from the original, raw data. Estimated travel times are based on data presented in Table 1. Obviously, more analysis needs to be performed on the composite materials. Additional data also needs to be taken on di€erent types of samples to attempt to quantify speci®c e€ects. 5. Conclusions and recommendations

Fig. 11. Reference signal (top) and through-transmitted signal (bottom) for Ti±Gr±Ti sample, 25 MHz (nominal) input center frequency.

A discrete-time H1 deconvolution ®lter based on the bounded real lemma has been presented. The application of the ®lter to nondestructive evaluation (NDE) of materials was discussed and simulation results for this class of applications were presented. The ®lter was applied to ultrasonic NDE data. Results obtained from the H1 deconvolution ®lter are compared to results obtained from a Wiener deconvolution ®lter. The H1 ®lter provides exact results in the absence of noise, as does the Wiener ®lter. In the presence of

Fig. 12. Deconvolved transmission sequence from Ti±Gr±Ti sample, 25 MHz center frequency input: (a) H1 ®lter (D=0.5, =10): (b) Wiener ®lter.

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

measurement noise, the H1 ®lter results were comparable to results obtained from the Wiener ®lter. In some of the test data analysis the H1 ®lter appeared to be slightly more sensitive than the Wiener ®lter, as it responded to events apparent in the raw data which were unnoticeable in the Wiener ®lter results. Overall, the performance of the H1 ®lter is surprisingly good in the presence of noise. The Wiener ®lter has access to all available data, for every estimate; it should provide excellent noise suppression. The H1 ®lter has access only to the estimate at the previous time step and the measured data at the current time step; this is a signi®cant handicap in the presence of random noise. The H1 deconvolution ®lter shows considerable promise in NDE applications. The primary e€ect which hinders its performance appears to be random measurement noise. It is recommended that the H1 ®lter design be modi®ed to allow smoothing (use of several future data points to estimate the system states and inputs at the current time). This should allow the ®lter to obtain considerably more accurate results in the presence of measurement noise. Hybrid approaches, incorporating the features of both L2 and H1 ®lters, should also be investigated in order to improve the deconvolution results in the presence of noise.

333

y…k ‡ 1† ˆ C…k ‡ 1†x…k ‡ 1† ‡ D…k ‡ 1†w…k ‡ 1† ˆ C…k ‡ 1†‰A…k†x…k† ‡ B1 …k†u…k† ‡B2 …k†w…k†Š ‡ D…k ‡ 1†w…k ‡ 1† the following error expressions are obtained e…k ‡ 1† ˆ A^ …k†e…k† h i ‡ A^ …k† ÿ A…k† ‡ B^ 2 …k†C…k ‡ 1†A…k† x…k† h i ‡ B^ 2 …k†C…k ‡ 1†B1 …k† ÿ B1 …k† uu …k† h i ‡ B^ 1 …k† ‡ B^ 2 …k†C…k ‡ 1†B1 …k† ÿ B1 …k† u0 …k† h i ‡ B^ 2 …k†C…k ‡ 1†B2 …k† ÿ B2 …k† B^ 2 …k†D…k ‡ 1†   w…k† w… k ‡ 1 † …A2†

"…k† ˆ C^ …k†e…k† h i ‡ C^ …k† ‡ D^ 2 …k†C…k ‡ 1†A…k† x…k† h i ‡ D^ 1 …k† ‡ D^ 2 …k†C…k ‡ 1†B1 …k† ÿ Im u0 …k† h i ‡ D^ 2 …k†C…k ‡ 1†B1 …k† ÿ Im uu …k† h i ‡ D^ 2 …k†C…k ‡ 1†B2 …k† D^ 2 …k†D…k ‡ 1†   w…k† w…k ‡ 1†

Acknowledgements

…A3†

The authors would like to thank Je€ Daniels for acquiring the data presented in this paper and the Boeing Aircraft Company for providing the GFRE samples.

In order to eliminate the ®lter errors' dependence on x(k) and u0(k), we choose

Appendix: H1 deconvolution ®lter derivation

A^ …k† C^ …k† D^ 1 …k† B^ 1 …k†

The state estimation error, e(k), and the input estimation error, "(k), are de®ned as follows:

The estimation error dynamics of the ®lter can now be written as

e…k† "…k†

ˆ x^ …k† ÿ x…k† ˆ u^ …k† ÿ u…k†

…A1†

Substituting Eqs. (5) and (6) into Eqs. (A1) and noting that

ˆ ˆ ˆ ˆ

A…k† ÿ B^ 2 …k†C…k ‡ 1†A…k† ÿD^ 2 …k†C…k ‡ 1†A…k† Im ÿ D^ 2 …k†C…k ‡ 1†B1 …k† B1 …k† ÿ B^ 2 …k†C…k ‡ 1†B1 …k†

e…k ‡ 1† ˆ Ae …k†e…k† ‡ Be …k†…k† "…k† ˆ Ce …k†e…k† ‡ De …k†…k† where

…A4†

…A5†

334

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335

2

3

uu …k† 7 w…k† 5 w…k ‡ 1† Ae …k† ˆ A…k† ÿ B^ 2 …k†C…k ‡ 1†A…k†

such that

6 …k† ˆ 4

Be …k† ˆ B^ 2 …k†DA …k† ÿ BA …k† Ce …k† ˆ ÿD^ 2 …k†C…k ‡ 1†A…k†

R2 …k† ˆ 2 I ÿ C…k†Q…k†CT …k† ÿ D…k†DT …k† > 0 …A10†

…A6†

Applying the Bounded Real Lemma to Eqs. (A5) and assuming that

De …k† ˆ D^ 2 …k†DA …k† ÿ L

R…k† ˆ DA …k†DTA …k†

with

‡ C…k ‡ 1†A…k†Q…k†AT …k†CT …k ‡ 1† > 0 …A11†

B…k† ˆ ‰B1 …k† B2 …k†Š BA …k† ˆ ‰B…k† 0Š

it can be determined that

DA …k† ˆ ‰C…k ‡ 1†B…k† D…k ‡ 1†Š L ˆ ‰ I m 0 0Š We now apply the bounded real lemma for timevarying systems [8] to Eqs. (A5), to minimize the error terms in an H1 sense. This lemma is stated as follows. A.1. Bounded real lemma Consider the following linear discrete time, time-varying system, x…k ‡ 1† ˆ A…k†x…k† ‡ B…k†u…k† y…k† ˆ C…k†x…k† ‡ D…k†u…k†; x…t0 † ˆ x0

…A7†

For a positive , the necessary and sucient conditions for k y k22 < 2 SUP 2 T 06ˆu2L2 ‰t0 ;tf Š k u k2 ‡x …t0 †S2 x…t0 †

R2 …k† ˆ 2 L ÿ De …k†DTe …k† ÿ Ce …k†Q…k†CTe …k† h i ˆ R1 …k† ÿ D^ 2 …k† ÿ BT …k†CT …k ‡ 1†Rÿ1 …k† h iT R…k† D^ 2 …k† ÿ BT …k†CT …k ‡ 1†Rÿ1 …k† …A12† where ÿ
R1 …k† ˆ 2 ÿ 1 L ‡ BT …k†CT …k ‡ 1†Rÿ1 …k† C…k ‡ 1†B…k†: …A13† The necessary condition for (A8) is that R2(k)>0. It is apparent from (A12) that this condition can be satis®ed by choosing D^ 2 …k† ˆ BT …k†CT …k ‡ 1†Rÿ1 …k†:

…A14†

…A8†

It can also be determined (following Ref. [9]) that

where S2 is a symmetric positive de®nite matrix, are that there exists a symmetric solution Q(k) for k2[t0,tfÿ1] to the Riccati di€erence equation

 Ae …k†Q…k†ATe …k† ‡ Be …k†BTe …k† ‡ Ae …k†Q…k†CTe …k†    ‡Be …k†DTe …k† R…k† Ce …k†Q…k†ATe …k† ‡ De …k†BTe …k† h i h iT ˆ J…k† ‡ B^ 2 …k† ÿ Ks …k† R…k† B^ 2 …k† ÿ Ks …k†

Q…k ‡ 1† ˆ A…k†Q…k†AT …k† ‡ B…k†BT …k† ÿ
‡ A…k†Q…k†CT …k† ‡ B…k†DT …k† Rÿ1 2 …k† ÿ
T T C…k†Q…k†A …k† ‡ D…k†B …k† Q…t0 † ˆ Sÿ1 2 …A9†

…A15† where  Ks …k† ˆ A…k†Q…k†AT …k†CT …k ‡ 1†  ‡ BA …k†DTA …k† Rÿ1 …k†

…A16†

T.C. Hanshaw et al. / ISA Transactions 38 (1999) 323±335 T J…k† ˆ A…k†Q…k†AT …k† ‡ BA …k†RTÿ1 1 BA …k†   T A…k†Q…k†AT …k†CT …k ‡ 1† ‡ BA …k†Tÿ1 1 DA … k †   T ÿ1 T Tÿ1 2 …k† C…k ‡ 1†A…k†Q…k†A …k† ‡ DA …k†T1 BA …k†

…A17†

T1 and T2(k) are as given in Eqs. (11). It can be seen that J(k) is equal to Q(k+1), from Eq. (10). We choose B^ 2 …k† ˆ Ks …k†. References [1] J. KrautkraÈmer, H. KrautkraÈmer, Ultrasonic testing of materials, 4th Edition, Springer Verlag, New York, 1990. [2] G. Hayward, J. Lewis, Comparison of some non-adaptive deconvolution techniques for resolution enhancement of ultrasonic data, Ultrasonics 27 (5) (1989) 155±164. [3] S.K. Sin, C.H. Chen, A comparison of deconvolution techniques for the ultrasonic nondestructive evaluation of materials, IEEE Transactions on Image Processing 1 (1) (1992) 3±10. [4] J.M. Mendel, Optimal seismic deconvolution Ð an estimation-based approach, Academic Press, New York, 1983.

335

[5] M.J. Grimble, H1 ®xed-lag smoothing ®lter for scalar systems, IEEE Transactions on Signal Processing 39 (9) (1991) 1955±1963. [6] K.M. Nagpal, P.P. Khargonekar, Filtering and smoothing in an H1 setting, IEEE Transactions on Automatic Control 36 (32) (1991) 152±166. [7] J.G. Doyle, B.A. Francis, A.R. Tannenbaum, Feedback control theory, Macmillan Publishing Company, New York, 1992. [8] X. Yu, C.S. Hsu, Reduced order H1 ®lter design for discrete time-varying systems, International Journal of Robust and Nonlinear Control 7 (1997) 797±809. [9] X. Yu, C.S. Hsu, R. Bamberger, S.J. Reeves, H1 deconvolution ®lter design and its application in image restoration, Proceedings of the International Conference on Acoustics, Speech and Signal Processing 4 (1995) 2611± 2614. [10] J.-N. Juang, R.S. Pappa, An eigensystem realization algorithm for modal parameter identi®cation and model reduction, Journal of Guidance, Control and Dynamics 8 (5) (1985) 620±627. [11] J.-N. Juang, R. Pappa, E€ect of noise on modal parameters identi®ed by the eigensystem realization algorithm, Journal of Guidance, Control and Dynamics 9 (3) (1986) 294±303. [12] J.-N. Juang, Applied system identi®cation, PTR PrenticeHall, Englewood Cli€s, New Jersey, 1994.