Particle filter based noise removal method for acoustic emission signals

Mechanical Systems and Signal Processing 28 (2012) 63–77

Contents lists available at SciVerse ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Particle filter based noise removal method for acoustic emission signals Changjiang Zhou, Yunfeng Zhang n Department of Civil and Environmental Engineering, University of Maryland, College Park, MD 20742, USA

a r t i c l e i n f o

abstract

Article history: Received 1 December 2010 Received in revised form 27 May 2011 Accepted 11 August 2011 Available online 8 September 2011

This paper discusses the application of a statistical noise removal technique – Rao– Blackwellised particle filter (RBPF) for signal to noise ratio (SNR) enhancement of acoustic emission (AE) signals. RBPF is a recursive Bayesian method for dynamic system state estimation. Compared with other signal filtering methods, RBPF offers the advantage of broad-band signal cleansing by directly modeling the internal dynamics of the concerned physical system and statistical characteristics of the signal noise. In doing so, signal filtering can be related to the dynamic characteristics of the underlying physical system, rather than a purely mathematical operation. RBPF also outperforms a few other statistical signal filtering methods such as Kalman filter, with the ability of handling nonlinear system and non-Gaussian noise problems. Another feature of RBPF is its ability to allow real-time on-board signal processing. In this paper, moment tensor analysis was performed first to generate simulated baseline AE signals. The simulated AE signal was subsequently superimposed with noise to demonstrate the effectiveness of the RBPF method in filtering noise in the AE signals. The results show that the performance of the RBPF in SNR enhancement is very promising. Finally, RBPF is also applied to real AE data obtained from experimental tests and apparent improvement to the SNR in AE feature study is observed. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Acoustic emission Filtering Particle filter Sensor Signal noise ratio Signal processing

1. Introduction Acoustic emission (AE) has been demonstrated as a useful nondestructive testing (NDT) technique for local damage detection in structural members made of various types of materials. AE is the elastic wave generated by sudden energy release within a material, which provides real-time information on damage progression within a structure. Since most of the AE sources are damage-related, AE monitoring can be effectively used to diagnose impending structural material failure. The frequency content and amplitude of the measured AE signal depend on the type of AE source, structural geometry, material properties as well as the frequency response characteristics of the adopted AE sensor. Different from ultrasonic testing, which excites elastic waves into a solid, AE sensor passively listens to the signals generated by crack initiation and progression within the monitored structure. The understanding of AE signals associated with different damage mechanisms is crucial to signal interpretation and decision making. In general, the AE signal analysis methods can be classified into two categories—the parameter based technique and waveform based technique. Instead of recording the AE waveform, the first method (also known as ring-down count test) only records AE features such as cumulative AE counts and RA value. Recently, with advancement in broad-band AE sensors, waveform based AE signal

n

Corresponding author. Tel.: þ1 301 405 1955. E-mail addresses: [email protected] (C. Zhou), [email protected] (Y. Zhang).

0888-3270/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2011.08.004

64

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

interpretation techniques have received growing attention with an aim of more accurately quantifying the crack source characteristics. Regardless of the method used, noise removal is essential to the quality of the information extracted from the measured AE signal. Broad-band AE data measurement is subject to noise contamination, which means other non-related events can be mixed up with the desired AE signal. Depending on noise characteristics, two major signal processing techniques are available to filter out the noise content from the original signal: frequency domain and time-space domain methods. If the noise and desired signal contents do not overlap with each other in the frequency domain, common approach is to filter out the noise through frequency contents selection, such as applying a low pass filter, or band-pass filter, etc. In the second approach, for example, white noise, which spreads over the entire frequency bandwidth, often cannot be removed in a traditional way through frequency bandwidth selection. Instead, a state space approach based on the statistical characteristics of the noise can usually be employed for noise removal. Examples of such a filter design in time domain are: Wiener filter [1], Kalman filter [2], Savitzky–Golay smoothing filter [3], etc. For the best performance of the Wiener filter, both signal and additive noise need to be a linear stationary process. If the process is non-stationary, Wiener filter cannot be used. Kalman filter is a very popular method for noise removal, which can be used for non-stationary signal process. However, to ensure accuracy when using the Kalman filter, the underlying dynamic system needs to be linear and the noise has to be a Gaussian-white noise process. Although Extended Kalman filter [4] has been proposed for nonlinear system application, the estimation problem still needs to be linearized about the predicted state and requires a Gaussiandistribution for the noise. In order to have more flexibility in dealing with nonlinear and non-Gaussian-distribution cases, a new signal filtering technique—particle filter has been recently developed [5]. Particle filter does not require the underlying physical system to be a linear system nor require Gaussian-distributed noise. It is a state space approach through a recursive Bayesian method for dynamic state estimation, in which the probability density function of the state vector is constructed based on available system information (e.g., a model describing the physical system). This paper deals with the application of Rao–Blackwellised particle filter (RBPF) for signal to noise ratio (SNR) enhancement of AE signal. Moment tensor analysis is performed first to generate simulated baseline AE signals. The simulated AE signal is subsequently superimposed with noise to demonstrate the effectiveness of the RBPF method in filtering noise in AE signals. The results show that the performance of RBPF in SNR enhancement is very promising for AE signal cleansing. Finally, RBPF is also applied to real AE data obtained in field and lab tests and apparent improvement in SNR is observed from a feature study of RA value and cumulative counts versus time of AE signals.

2. Overview of particle filter 2.1. Particle filter Particle filter, also referred to as sequential Monte Carlo method, is a simulation and statistical model-based estimation method [4]. It can be used to construct the posterior probability distribution of the signal based on available measurements and the assumption of the noise probability distribution. Particle filters are similar to Markov chain Monte Carlo (MCMC) batch methods and only estimate the current state based on the available observations (samples). However, most of the MCMC methods would require a fixed batch of dataset to model the posterior density function, which is not suitable for online filtering [6]. In order to implement the particle filter, a state space model is needed Xk ¼ fk1 ðXk1 ,Uk1 Þ

ð1-aÞ

Zk ¼ hk ðXk ,Vk Þ

ð1-bÞ

where Xk, Zk are the state vector and the measurement vector of the dynamic system, respectively, at time step k. f and h are state function and observation function, respectively, while Uk, Vk are noise process for the state and observation, respectively. Capital letter here denotes the stochastic process while lowercase letter denotes the value of the stochastic variables, for example, zk is the value of Zk at time step k. Particle filter aims to construct the posterior probability distribution pðXk 9z1:k Þ based only on the observed data. According to the Chapman–Kolmogorov equation Z pðXk 9z1:k1 Þ ¼ pðXk 9Xk1 ÞpðXk1 9z1:k1 ÞdXk1 ð2Þ where pðXk 9Xk1 Þ is the transition probability given by the state space. According to the recursive Bayesian estimation pðXk 9z1:k Þ ¼

pðzk 9Xk ÞpðXk 9z1:k1 Þ pðzk 9z1:k1 Þ

ð3Þ

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

where the denominator Z pðzk 9z1:k1 Þ ¼ pðzk 9Xk ÞpðXk 9z1:k1 ÞdXk

65

ð4Þ

One can first sample X0i  pðX0 Þ, i ¼ 1 : NS , where NS is the number of particles (differently weighted samples of the i distribution) to be employed, then for each particle X0i , sample X~ 1  pðX1 9X0i Þ. For each particle, based on the measurement equation expressed above as Eq. (1-b), starting at k¼1, the importance weights can be measured with i ~ i ¼ pðZk 9X~ k Þ w

ð5Þ

The next step is to normalize the particle weights obtained in the previous step, so that the summation of the total weights of the particles equals one and the weight associated with each particle can represent the probability distribution ~i w wi ¼ PN S

i¼1

ð6Þ

~i w

i Theses weights are associated with X~ k . i Next resample NS particles from ðX~ k ,i ¼ 1,. . .,NS Þ to obtain ðXki ,i ¼ 1,. . .,NS Þ. Based on the probability density distribution of Xki , the optimal estimation of the state X at time step k can be obtained as follows:

X^ k ¼ Xki wi

ð7Þ

The above process is then repeated for k ¼2, 3, 4, etc. 2.2. Rao–Blackwellised particle filter For AE signals, which are normally recorded at high sampling rates (on the order of mega samples per second), particle filter can become cumbersome due to the need for large samples to construct the posterior probability distribution. This is especially severe in high dimensional state-spaces and in real-time monitoring application. In order to more efficiently implement the particle filters in AE signal monitoring, it is thus important to marginalize out some variables so that the size of the state space variable vector can be reduced. By partitioning the state-space Xk into two sub-spaces, Rk and Sk , based on conditional probability, it yields pðR1:k 9S1:k ,Z1:k ÞpðS1:k ,Z1:k Þ pðR1:k ,S1:k ,Z1:k Þ pðS1:k ,Z1:k Þ ¼ ¼ pðR1:k 9S1:k ,Z1:k Þ pðZ1:k Þ pðZ1:k Þ pðZ1:k Þ ¼ pðR1:k 9S1:k ,Z1:k ÞpðS1:k 9Z1:k Þ

pðR1:k ,S1:k 9Z1:k Þ ¼

ð8Þ

If pðR1:k 9S1:k ,Z1:k Þ can be analytically updated, the only part that requires sufficient number of particles would be pðS1:k 9Z1:k Þ. As the state space S1:k is a sub-space of the original state space, fewer particles are needed to implement the particle filter. This theorem is referred to as the Rao–Blackwellised particle filter (RBPF) [7]. For particle filters including the RBPF, the restrictions on noise characteristics are not as stringent as the Kalman filter, namely, they do not have to be Gaussian-distributed, nor the state space needs to be linear. In the following discussions for the concerned AE signal cleansing application, noise with Gaussian-distribution characteristics is assumed as examples to illustrate the methodology. For simplicity, the state space system is also assumed to be linear, while the nonlinear feature will be accounted for by the observation equation and taken care of by the particles in the way of a jump Markov linear system, which will be explained in detail later in the article. 3. System model for particle filter based signal cleansing 3.1. Model formulation in state space Earthquakes are caused by the sudden release of energy during the fracture of stressed rock within the Earth’s crust. This phenomenon is similar to that process, which occurs in materials under load, and although they take place on very different scales, these two phenomena – earthquakes in geophysics and damage in structural materials – have similarities. In both cases, there is a release of elastic energy from sources located inside a medium [15]. In seismology field, an amplitude modulated cyclic waveform is usually used to represent the seismic wavelet [11]. Naturally, because of the similarity between seismic signals and acoustic emission signals, this seismic wavelet model can also be considered as a good candidate for representing the AE events in structural field. For tZt0, this model [8] is expressed as Aðtt0 Þ ¼ A1 ðtt0 Þsin½oðtt0 Þ

ð9Þ

where A1 ðtt0 Þ is the amplitude of the AE source wavelet, and t0 is the wavelet arrival time; o is the dominant angular frequency. In a similar way to those used by other researchers [9–12] in processing the seismology data, a stochastic process X2 ðtÞ is adopted to represent the seismic wavelet’s amplitude while another stochastic process Xw ðtÞ is used to represent the seismic wavelet. The process can be rewritten into the following form: Xw ðtÞ ¼ X2 ðtÞsin½oðtÞtdðtÞ

ð10Þ

66

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

where dðtÞ is the phase angle. As the Gauss–Markov process can be used to model many physical phenomena, it is also chosen here to model the seismic wavelet amplitude X2 ðtÞ, which is expressed as follows in a discrete form: ð11Þ X2 ðk þ 1Þ ¼ aX2 X2 ðkÞ þ bX2 w2 ðkÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi bD 2 bD , bX2 ¼ s ð1e Þ, and w2 is a zero-mean, timewise-uncorrelated, unit-variance sequence with Gaussian where aX2 ¼ e probability distribution function; b ¼ 1=Tc , where Tc is the time constant of the Gauss–Markov process; s is the standard deviation for the process;D is the sampling time interval. Assuming that the noise part of the state space can also be modeled as a Gauss–Markov process, it yields X1 ðk þ 1Þ ¼ aX1 X1 ðkÞ þ bX1 w1 ðkÞ

ð12Þ

where aX1 , bX1 are the Gauss–Markov process parameters, and w1 is a zero-mean, timewise-uncorrelated, unit-variance sequence with Gaussian probability distribution function independent of w2 . All of the above parameters can be determined from the noise part of the AE signal. Details can be seen in the RBPF algorithm outlined in Section 3.2. Using the above-described model, the state space system is formulated as: #" #" " # " # " # a X1 0 bX1 0 X1 ðkþ 1Þ X1 ðkÞ w1 ðkÞ X¼ ¼ þ ð13Þ 0 aX2 0 bX2 X2 ðkþ 1Þ X2 ðkÞ w2 ðkÞ Accordingly, the measurement function can be constructed as: " # X1 ðkÞ ZðkÞ ¼ ½1 sinðok kDt þ dk Þ þ vk X2 ðkÞ

ð14Þ

where ok is the dominant frequency at time step k, and vk is the measurement noise. Using a matrix form to represent the equations yields Xk þ 1 ¼ Fk Xk þ Bk wk

ð15-aÞ

Zk ¼ Hk Xk þ vk

ð15-bÞ

The above state function and measurement function will be used in steps 2, 3, 4 in the RBPF Algorithm described in Section 3.2. It should be noted that the denoised signal is the X2 in the state vector. The reconstructed signal would be the product of the second component of the observation function (y~ k , estimated by particle filter, see Algorithm 2: Particle filter updating for details) and X2. However, since there is no fixed dominant frequency, a set of particles of frequency oik with higher likelihood will be used in order to reduce the inaccuracy caused by a specifically selected frequency. In this study, the dynamic and measurement functions take the following form for each particle: i Xki ¼ Fk1 ðyik ÞXk1 þ Gk1 ðyik Þuik1

ð16-aÞ

Zki ¼ Hk ðyik ÞXki þ vik

ð16-bÞ

It is seen now that for each particle the state and the measurement functions are linear. In other words, in the proposed RBPF, standard Kalman filter updating is used for state updating while its non-linearity capability is realized through particle filter in the measurement function. As the RBPF model is extended from seismology field to AE signals in engineering applications, the frequency range of interest is substantially increased from a few Hertz to several kilo Hertz [16]. Therefore, if the same grid based particle filter is still used for system formulation as in the seismology area, an overly large number of particles oik would be required to obtain the same level of accuracy over the entire frequency range. In order to reduce the total number of particles, the particles (frequency grids) are redistributed into a logarithm form in this study rather than a uniform distribution adopted in the seismology study by others [9–12]. This implies that logðoik Þ follows a uniform distribution over the desired frequency range. The general procedure for implementing particle filter in AE signal noise removal can be summarized as follows: a) Specify the state function (Eq. (1-a)) and the measurement function (Eq. (1-b)). b) Choose a proper stochastic process to represent the state vector and the measurement vector as well as the noise characteristics. c) Initialize the particle filter and then perform signal filtering following Algorithm 1 described next.

3.2. Pseudo code To illustrate the algorithm of the RBPF, the procedure is presented in the following pseudo code:

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

67

Algorithm 1. RBPF updating NS ~ k ,Z j gNS  ½fXkj ,wjk gj ¼ 1  ¼ RBPF½fH k1 j ¼ 1

 Step 1: Initialize the Kalman filter by letting X1 ð0Þ ¼ Z0 , X2 ð0Þ ¼ 0; set error covariance P090 ¼ ½R 0; 0 R, where R is



the initial value of the error covariance, If the exact initial value of the state vector is known, R can be set as zero, otherwise, R can be set to be a relatively large value. The writers’ experience is that selection of a relatively large R value (e.g., 50 times the signal amplitude) generally yields a satisfactory result. From the writers’ research, it is found that R value selection has little effect on the filtering results once the propagation converges. For Time_i¼ 1, Step 2: Propagate the state estimate extrapolation and error covariance matrix one step ahead: X~ k9k1 ¼ Fk1 X~ k19k1 P~ k9k1 ¼ Fk1 Pk19k1 FTk1 þQk1 Here, Fk can be initialized based on prior knowledge of the state function. It can be time independent (as used in this paper) or time dependent, which needs updating or propagating at each time step. In this paper, Fk is set to be [a1 0; 0 a2] for k¼1, where ai , i ¼ 1,2, are the parameters in the Gauss–Markov process as shown in Eqs. (11) and (12) (different values for different state variables, i.e. a1, a2). The appropriate value for a1 and a2 can be obtained from the signal sequence (with noise) and the noise part (which can be collected by taking samples without the presence of the desired signals). Qk is the covariance of the state noise process: uk ¼ Bk wk  ð0,Qk Þ, in other words,

T w Bk

Qk ¼ B k s

¼

b1

0

0

b2

!

s21

0

0

s22

!

b1

0

0

b2

!T ¼

s21 b21

0

0

s22 b22

!

It is assumed here that the variance of noise does not change over time (i.e., time invariant random process). Thus, Qk is the same as Q, which could be obtained from statistical analysis of the noise part of the signal or based on prior knowledge of the noise process. In this paper, the Q value (matrix) is initialized based on the noise part of the signal. The general procedure to determine Q value employed in this study is given as follows: (a) Measure the noise portion of the signal without the presence of the desired signal (i.e., collecting ambient noise signal). (b) Calculate the autocorrelation of the sampled noise signal. (c) Perform a curve fitting of the autocorrelation data using the model, s2 eb9t9 [11], in which s and b (1/Tc, the time constant) are going to be determined from regression. qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi After obtaining sX1 and bX1 , calculate aX1 ¼ ebX1 D and bX1 ¼ sX1 ð1e2bX1 D Þ. It is seen that for the noise process, when b is very large (implying that signal is less correlated), bX1  sX1 , bX1 ðb1 Þ can be initialized as the value of standard deviation of the noise signal [11]. For b2, some trial and error study has to be done to obtain an acceptable value. It was found in this study that setting the b2 value equal to bX1 ðb1 Þ offers a good starting point. N i  Step 3: Particle update: ½fH~ k ,wjk19k1 gN  ¼ PF½fwik19k1 ,Zk1 g *. S

S

i¼1

j¼1

 Step 4: Use the updated measurement matrix to recalculate measurement extrapolation ~ k X~ Z~ k ¼ H k9k1 and variance of innovation ~T ~ kP S~ k ¼ H k9k1 H k þ Rk



where Rk is the variance of the observation noise in the measurement function, vk, i.e. vk  ð0,Rk Þ. Technically speaking, this Rk can be updated from measured data by calculating the variance of the observation noise or obtained from sensor noise specification. In this paper, the variance of the observation noise is assumed to be a constant and set to a relatively small value. Step 5: Innovation can be updated as: I~ k ¼ Zk Z~ k T 1

 Step 6: Kalman gain matrix update: K~ k ¼ Pk9k1 H~ k S~ k .  Step 7: State estimation update: X~ k9k ¼ X~ k9k1 þ K~ k I~k .

68

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

 Step 8: Error covariance update: P~ k9k ¼ ½IK~ k I~k Pk9k1 where I is the identity matrix in the above equation.  Step 9: Repeat the process from Time_i ¼ 2 : k (up to date).

*where the PF function can be implemented as the following procedure: Algorithm 2. Particle filter updating ~ k ,wj g ½fH k19k1

NS j¼1

i  ¼ PF½fwik19k1 ,Zk1 g

NS i¼1



 Sample oik  pðok Þ, i ¼ 1 : NS  For yik ¼ sinðoik kDÞ, calculate the measurement extrapolation i Z~ k ¼ Hki X~ k9k1

and the variance of the innovation Sik ¼ Hki Pk9k1 ðHki ÞT þ Rk

 Calculate the predictive density [13] i

pðZk 9yik Þ ¼ pðZk 9Z1:k1 ,yi1:k Þ ¼ NðZk ; Z~ k ,Sik Þ

 Calculate the importance weight of each particle wik9k1 ¼

NS X

wjk19k1 pðyik 9yjk1 Þ

j¼1

wik9k1 pðZk 9yik Þ wik9k ¼ PN j j S j ¼ 1 wk9k1 pðZk 9yk Þ

 Obtain estimate of y~ k by its mean value y~ k ¼ E½yk 9z1:k  ¼

NS X

wjk9k yjk

j¼1

 Update measurement matrix ~ k ¼ ½1 H

y~ k :

4. RBPF based signal cleansing: benefits and examples 4.1. Advantages of RBPF in AE signal cleansing 4.1.1. Real-time processing feature The RBPF, which is a state space time domain filtering technique, is updated at every time step when a new data point is received. It does not depend on the previous data except for the very last one. Therefore, for execution, there is no special requirement for data record length and sampling frequency. For example, Fig. 1 shows a measured AE data collected from a pencil lead break test, which is adopted as a common approach to AE event simulation. As seen in Fig. 1, the filtering effect is not affected by the length of the received data at the time of processing. Future data is thus not required for signal cleansing with RBPF. This feature is not available from conventional filters, which may require the total dataset for processing. For onboard processing (e.g., signal cleansing on wireless sensor node), this is a very appealing feature, as it enables real-time data processing on the fly, without worrying about how to choose a proper length for dataset to be processed. 4.1.2. Numerical simulation example 4.1.2.1. Moment tensor theory. Many AE signal sources can be analyzed with moment tensor theory, especially for the cracking type sources [17]. Good agreement has been found between experimental data and theoretical synthetic data. In

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

69

Original data (100%)

Signal Amplitude

20% data length (filtered) 50% data length (filtered) 70% data length (filtered) 100% data length (filtered)

0

0.05

0.1 0.15 Time (ms)

0.2

0.25

Fig. 1. Demonstration of the real-time processing feature of RBPF based signal filtering.

this study, AE signals generated from moment sensor analysis will be used to represent AE signals when performing parametric study of AE signal cleansing with RBPF. Usually, moment tensor analysis employs a model combining multiple dipoles and double couples. Without loss of generality, it is assumed that crack orientation is parallel to one of the coordinate axes of the coordinate system. As in the generalized AE theory [14], AE waves due to cracking can be expressed as: uk ðx,tÞ ¼ Gkp,q ðx,y,tÞSðtÞCpqij nj li DV

ð17Þ

where Cpqij are elastic constants, DV is the crack volume, n is the crack normal and l is the crack motion vector; SðtÞ is the source mechanism; Gkp,q is the derivative of Green’s function Gkp with respect to coordinate q. Denoting Cpqij nj li DV as Mpq , also referred to as moment tensor, yields uk ðx,tÞ ¼ Gkp,q ðx,y,tÞSðtÞMpq

ð18Þ

With different combinations of moment tensors and Green’s functions, different types of cracks can be modeled with Eq. (18). For illustration purpose, a case of shear type crack is considered, which represents a typical cracking mode—fatigue in cyclically loaded structures. Also, instead of displacement, strain signal is going to be synthesized here since experimental AE data in this study are all collected from stress wave AE sensors. 4.1.2.2. Synthetic condition. Assume a crack happens in this way: crack surfaces F þ and F  are on the x1–x2 plane, normal to the x3-axis; the surface F þ moves in horizontal and perpendicular to the normal vector n when the crack opens. If the normal vector is defined as n¼ (0, 0, 1), the motion vector would be l ¼(1, 0, 0). The source mechanism is assumed as follows: S0 ðtÞ ¼ sin4 ðpt=TrÞ, S0 ðtÞ ¼ 0,

t 4 Tr

t rTr

ð19-aÞ ð19-bÞ

where Tr¼ 2 ms is the rise time of this cracking event. Providing there is a subsurface crack with a depth of 0.002 m, and sensor for AE signal measurement is located 10 cm (4 in.) away from the crack projection on the surface, synthetic AE signal obtained from moment tensor analysis is shown in Fig. 2(a). The sampling frequency is set to be 10 mega samples per second and the total sampling duration is 102.4 ms. To demonstrate the effectiveness of RBPF in signal noise removal, a noise process with a Gaussian-distribution probability density is added to the simulated AE signal as shown in Fig. 2(b). The additive noise process has a variance of 0.0012 and a time constant of 2e6 ms. Fig. 2(b) shows the signal with the additive noise. It is seen from Fig. 2(b) that the SN ratio is quite low and removal of the noise is highly desired for any meaningful post-processing of AE data such as counts rate calculation. Fig. 3(a) shows the filtered signal after applying a conventional band-pass filter with a pass band between 0.5 kHz and 2 MHz. The filter used is a second order Butterworth filter. For comparison purpose, Fig. 3(b) shows the signal filtered with the RBPF method. The frequency particles are specified by 500 points with logarithm distributed across the specified frequency range. From Fig. 3(b), it is seen that the AE signal filtered with RBPF is much cleaner and more closely resembles the original signal (i.e., the one without noise in Fig. 2(a)) than the one filtered with the Butterworth filter. This means that in time domain, RBPF outperforms the band-pass filter. For AE signal, its frequency content may disclose very important information about the AE source. Therefore, Fig. 4(a)–(d) plots the frequency (amplitude) spectrum of the AE signals with or without filtering. It is again observed that filtering with RBPF gives better results than the Butterworth filter. Furthermore, in the pass band of the Butterworth filter, the frequency spectrum quality of the RBPF-filtered signal is better than the one filtered with the Butterworth filter. Clearly, a band-pass filter is unable to filter out a noise content, which

70

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

Amplitude (Volt)

2

1

0

−1 0

20

40

60

80

100

120

80

100

120

Time (μs)

Amplitude (Volt)

2

1

0

−1 0

20

40

60 Time (μs)

Fig. 2. Synthetic AE signal from moment tensor analysis: (a) original signal and (b) signal with additive noise.

Amplitude (Volt)

1.5 1 0.5 0 −0.5 −1 0

20

40

60

80

100

120

80

100

120

Time (μs)

Amplitude (Volt)

1.5 1 0.5 0 −0.5 −1 0

20

40

60 Time (μs)

Fig. 3. (a) Synthetic AE signal filtered with conventional filter (2nd-order Butterworth filter with a pass band between 0.5 KHz and 2 MHz) and (b) synthetic AE signal filtered with RBPF.

falls within its pass band frequency range. From this example, it is concluded that the RBPF can achieve a better filtering effect than the conventional filter when the noise content overlaps with the desired signal in certain frequency range. 4.2. Examples with measured data The following data were obtained from field tests on two steel highway bridges in South Korea. Due to the high noise level of the AE sensors observed in the field tests, SNR of the measured AE data was very poor. In order to extract useful

0.04

0.04

0.03

0.03 Amplitude

Amplitude

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

0.02

0.01

0.02

0.01

0 104

105 106 Frequency (Hz)

0.04

0.04

0.03

0.03 Amplitude

Amplitude

0 104

0.02

105 106 Frequency (Hz)

0.02

0.01

0.01

0 104

71

0 105

106

Frequency (Hz)

104

105

106

Frequency (Hz)

Fig. 4. Frequency spectra of synthetic AE signals: (a) original synthetic AE signal, (b) synthetic AE signal with additive noise, (c) filtered AE signal with additive noise (filtered with 2nd-order Butterworth band-pass filter) and (d) filtered AE signal with additive noise (filtered with RBPF).

information from the noisy AE signal sequence, the above-described RBPF method was applied to enhance the SNR of the measured data. The RBPF is specified by a frequency window of 0–500 KHz. The number of total particles was chosen to be 500; the noise is modeled as the Gauss–Markov process and is calculated to have a time constant of 0.00032 ms and a variance of 0.026 V2. The variance of the source wavelet is obtained from the signal and set to be 0.047 V2 and its time constant is 0.0344 ms. The sampling frequency is 2.5 MHz. Here the following assumptions are made: the initial wavelet amplitude is zero, x2 ð1Þ ¼ 0, and the first available measurement is set to be noise, i.e., x1 ð1Þ ¼ z1 .

4.2.1. Impact-induced AE signals In this field test, a total of six piezo paint AE sensors were installed on a steel box-girder bridge located near Kimpo, South Korea. Fig. 5(a) shows the signals measured by the piezo paint AE sensor mounted on the middle height of the web of the steel box girder. The signal was the response of the piezo paint AE sensor to an impact made with an impact hammer approximately 150 mm away from the sensor. It is seen in Fig. 5(a) that noise with a large amplitude dominated the measured signal. Further analysis of the measured data using FFT revealed that the noise was caused by the radio wave interference from a nearby communication tower. Fig. 5(b) shows the signal after implementing RBPF. It is seen from Fig. 5(b) that the signal after filtering is much cleaner than the original one. The SN ratio is greatly enhanced by the RBPF filtering. In order to ensure that the RBPF filtering procedure would not lead to loss of key information especially in the frequency range of interest, the frequency spectrum of the signal is shown in Fig. 6. It can be seen that the filtered signal still keeps the main frequency contents of the original signal while removing most of the broad band noise contents. Another available data was also collected from the same bridge test with different AE source condition: the simulated AE source was located 35 cm away from the sensor instead of 15 cm for the previous data. Fig. 7 shows a comparison of the original signal and the one after implementing the RBPF filtering. Again, RBPF did a fairly good job in cleansing the noisy AE

72

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

Amplitude (Volt)

1.5 1 0.5 0 −0.5 −1 0

0.5

1

1.5

2 Time (ms)

2.5

3

3.5

4

0

0.5

1

1.5

2 Time (ms)

2.5

3

3.5

4

Amplitude (Volt)

1.5 1 0.5 0 −0.5 −1

0.12

0.12

0.1

0.1

0.08

0.08 Amplitude

Amplitude

Fig. 5. Impact-induced AE signal (case I: 15 cm distance): (a) original signal with noise and (b) filtered signal with RBPF.

0.06

0.06

0.04

0.04

0.02

0.02

0 102

104

106

Frequency (Hz)

108

0 102

104

106

108

Frequency (Hz)

Fig. 6. Frequency spectra of impact-induced AE signal (case I: 15 cm distance): (a) original signal with noise and (b) filtered signal with RBPF.

data. Viewing the frequency spectrum shown in Fig. 8, the same conclusion made for Fig. 6 can be drawn for high-fidelity retaining of frequency contents of interest.

4.2.2. Vehicle-induced AE signal This field test was conducted on a steel I-girder bridge near Icheon, South Korea. Vehicles running on bridges are common AE sources in civil engineering. This type of AE signal is very useful because it is related to the bridge motion and itself is an indicator of the deformation experienced by bridge element subject to dynamic loading caused by moving vehicles. Therefore, obtaining a clean AE signal in a field test of bridge is of great interest to engineers. However, moving vehicles may also induce noise into measured AE signal, for instance, due to friction between tire and deck surface. These noises may spread over a frequency range of interest for diagnosing bridge condition using AE signal analysis. Therefore, the RBPF was employed again to test its performance in cleansing such data. Fig. 9 shows a signal sequence obtained when a truck with 25 ton gross vehicle weight (GVW) crossed the bridge at a speed of 80 km/h. The frequency domain performance is shown in Fig. 10. The RBPF is found to perform well in filtering data of this type.

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

73

Amplitude (Volt)

1 0.5 0 −0.5 −1 0

0.5

1

1.5

2

2.5

3

3.5

4

2.5

3

3.5

4

Time (ms)

Amplitude (Volt)

1 0.5 0 −0.5 −1 0

0.5

1

1.5

2 Time (ms)

0.15

0.15

0.10

0.10

Amplitude

Amplitude

Fig. 7. Impact-induced AE signal (case II: 35 cm distance): (a) original signal with noise and (b) filtered signal with RBPF.

0.05

0 102

0.05

0 104 106 Frequency (Hz)

108

102

104 106 Frequency (Hz)

108

Fig. 8. Frequency spectra of impact-induced AE signal (case II: 35 cm distance): (a) original signal with noise and (b) filtered signal with RBPF.

4.2.3. AE signal from fatigue test This fatigue test was conducted on a steel bridge orthotropic deck under cyclic fatigue loading. During the fatigue test, a piezo paint AE sensor was mounted on the specimen surface to monitor AE signals released by fatigue crack growth. Fig. 11(a) shows AE data collected during the fatigue test. It should be noted that the piezo paint sensor is close to an existing visible fatigue crack. Fig. 11(b) is the signal after denoising by RBPF. It can be seen that the RBPF performs well in filtering real AE data of this type. 4.3. Filter effect assessment The above examples show that AE signals after RBPF filtering have much improved SNR ratio. This is very beneficial to the AE signal feature study, which is commonly used in the AE field. The following section involves an AE feature study on the reconstructed signals with RBPF. In particular, results of RA value and cumulative AE counts versus time are presented next. 4.3.1. RA value RA value is an important AE feature in waveform based analysis, which is also known as the reciprocal of gradient in AE waveforms [18]. It is defined as the rise time (the time between the signal first exceeding the threshold value and the peak

74

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

Amplitude (Volt)

0.04 0.02 0 −0.02 −0.04 −5

−4

−3

−2

−1

0 1 Time (ms)

2

3

4

5

Amplitude (Volt)

0.04 0.02 0 −0.02 −0.04 −5

0

5

Time (ms) Fig. 9. Vehicle-induced AE signal measured on a steel I-girder bridge (vehicle speed ¼ 80 km/h): (a) original signal with noise and (b) filtered signal with RBPF.

x 10−3

6

5

5

4

4

Amplitude

Amplitude

6

3

3

2

2

1

1

0 102

104 Frequency (Hz)

106

x 10−3

0 102

104 Frequency (Hz)

106

Fig. 10. Frequency spectra of vehicle-induced AE signal measured on a steel I-girder bridge (vehicle speed ¼ 80 km/h): (a) original signal with noise and (b) filtered signal with RBPF.

amplitude arrival) divided by the peak amplitude of the signal [19]. It is believed that RA value is a very good index for classifying the fracture type: tensile type crack or shear type crack (or other types of crack); thus it is useful to identify failure modes of AE events. For example, it is found that lower RA values are more related to tensile type crack [18,20]. Since the baseline signal is known in the simulation part from moment tensor analysis, this synthetic AE data will be used here to study the effect of RBPF filtering on AE feature quality. In Fig. 12, the comparisons are based on such synthetic data. It is seen that the RA value (the reciprocal of the red dashed line slope) of the AE signal after applying the RBPF is much closer to the true value at a noise level of 50 mV (shown in Fig. 12(a)) if considering the RA value of the main portion of the signal. It is concluded that the RA value can be obtained more precisely with RBPF processed AE signals. 4.3.2. Cumulative AE counts Another example of examining the denoising effect of RBPF is presented by comparing the cumulative AE counts versus time. It is shown in Fig. 13 that the cumulative AE counts of a combined AE sequences (comprised of ten concatenated data sets each acquired from the different stages of a single fatigue test) at different threshold level. When the noise level is set

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

75

Amplitude (Volt)

0.15 0.10 0.05 0 0

0.05

0.1

0.15

0.2

0.25

0.15 Time (ms)

0.2

0.25

Time (ms)

Amplitude (Volt)

0.15 0.1 0.05 0 0

0.05

0.1

Amplitude (Volt)

Amplitude (Volt)

Fig. 11. AE signal measured from fatigue specimen: (a) original signal with noise and (b) filtered signal after denoising with RBPF.

1 0 −1 10

20

30 Time (μs)

40

50

1 0 −1 10

15

20

25

30

35

40

45

50

35

40

45

50

Amplitude (Volt)

Time (μs)

1 0 −1 10

15

20

25

30

Time (μs) Fig. 12. RA values of AE signals at a noise level of 50 mV: (a) original synthetic signal without noise, (b) filtered synthetic signal with additive noise (2nd-order Butterworth filtering with a pass band between 0.5 KHz and 2 MHz) and (c) filtered synthetic AE signal with additive noise (RBPF filtering). (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

to 10 mV, one can see that the noise spikes create a lot of counts in the analysis, which is not apparently true. In contrast, signals denoised with the RBPF algorithm exhibit much less variations and thus less error in cumulative AE counts calculation. It is concluded that the cumulative counts analysis of AE signals after using RBPF denoising seems to be more accurate due to less number of noisy spikes. For AE signals, their statistical values are more sensitive to spikes in post-processing since these spikes contaminate the AE signal, and pose significant difficulties in AE feature study, such as the RA value and cumulative counts analysis.

76

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

250 10mV(F) 10mV(O) 20mV(F) 20mV(O) 30mV(F) 30mV(O)

Cumulative Hits (log Σ N)

200

150

100

50

0 0

0.5

1

1.5

2

2.5

Time (ms)

Fig. 13. AE counts of signals versus time before (O) and after RBPF filtering (F).

The RBPF based denoising method is shown to be very promising in removing the noisy spikes and improving the SN ratio of AE signals. 5. Conclusions This paper presents a time-space filter known as Rao–Blackwellised particle filter (RBPF) and its application to denoising of acoustic emission (AE) signals. The RBPF can be used for both linear and nonlinear system model. Through properly choosing the model parameters for RBPF, it is found that the signal to noise (SN) ratio of AE signal can be considerably enhanced after implementing the RBPF based signal cleansing procedure. A validation test was conducted using simulated AE signals from moment tensor analysis. It is also found that owing to its time space feature, RBPF has the advantage of removing broad-band noise from the signal without corrupting the frequency content of interest in the original data. Real data obtained from field test on steel bridges and steel orthortropic deck specimen under fatigue loading was also used to demonstrate the effectiveness of RBPF in noise removal for AE data. The AE data include several common AE signal types in civil engineering. The results show that a substantial enhancement of SN ratio after filtering can be achieved using RBPF while the desirable frequency contents are retained. It is worth noting that the model utilized in this study was adapted from that originally derived for seismic data processing by other researchers. Considering the similarity and scale difference between seismic data and AE data from structural materials, this work extends the application of RBPF to ultrasonic AE signals through modification in model parameters, and subsequently verifies its effectiveness in AE signal cleansing. However, further research is still ongoing to design the optimal system model and noise characteristic process for AE signal cleansing.

Acknowledgments The work described in this paper was partially supported through a grant from the US National Science Foundation under Award nos. CMMI-0913857 and CMMI-0829327 (program manager: Dr. S.C. Liu). Partial financial support received from the University of Maryland and Minta Martin Aerospace Fund to conduct this research is also acknowledged. However, the opinions and conclusions expressed in this paper are solely those of the writers and do not necessarily reflect the views of the sponsors. References [1] N. Wiener, The Interpolation, Extrapolation and Smoothing of Stationary Time Series, Report of the Services 19, Research Project DIC-6037 MIT, February 1942. [2] R.E. Kalman, A new approach to linear filtering and prediction problems, Transactions of the ASME—Journal of Basic Engineering 82 (1960) 35–45. [3] A. Savitzky, M.J.E. Golay, Smoothing and differentiation of data by simplified least squares procedures, Analytical Chemistry 36 (8) (1964) 1627–1639. [4] M.S. Arulampalam, S. Maskell, N.J. Gordon, T. Clapp, A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking, IEEE Transactions on Signal Processing 50 (2) (2002) 174–189. [5] A. Doucet, N. Freitas, N.J. Gordon, An introduction to sequential Monte Carlo methods, in: A. Doucet, N. Freitas, N. Gordon (Eds.), Sequential Monte Carlo Methods in Practice (Statistics for Engineering and Information Science), Springer, 2001, pp. 3–14. [6] A. Doucet, N.J. Gordon, V. Krishnamurthy, Particle filters for state estimation of jump Markov linear systems, IEEE Transactions on Signal Processing 49 (3) (2001) 613–624.

C. Zhou, Y. Zhang / Mechanical Systems and Signal Processing 28 (2012) 63–77

77

[7] A. Doucet, N. Freitas, K. Murphy, S. Russell, Rao–Blackwellised particle filtering for dynamic Bayesian networks, in: C. Boutilier, M. Godszmidt (Eds.), Uncertainty in Artificial Intelligence, Morgan Kaufmann Publishers, 2000, pp. 176–183. [8] R.E. Sheriff, L.P. Geldart, Exploration Seismology, vol. 1, second ed., Cambridge University Press, Cambridge, UK, 1982, pp. 55. [9] E. Baziw, I. Weir-Jones, Application of Kalman filtering techniques for microseismic event detection, Pure and Applied Geophysics 159 (2002) 449–473. [10] E. Baziw, B. Nedilko, I. Weir-Jones, Microseismic event detection Kalman filter: derivation of noise covariance matrix and automated first break determination for accurate source location estimation, Pure and Applied Geophysics 161 (2004) 303–329. [11] E. Baziw, T.J. Ulrych, A Rao–Blackwellised type algorithm for passive seismic event detection, in: Proceedings of the Eighth Technical Meeting of the Consortium for the Development of Specialized Seismic Techniques (CDSST-2004), Vancouver, 14–15 December, Department of Geophysics, University of British Columbia, Vancouver, Canada, pp. 135–164. [12] E. Baziw, Real-time seismic signal enhancement utilizing a hybrid Rao–Blackwellised particle filter and hidden Markov model filter, IEEE Geoscience and Remote Sensing Letters 2 (4) (2005) 418–422. [13] N. Freita, Rao–Blackwellised particle filters for fault diagnosis, IEEE Aerospace Conference Proceedings 4 (2002) 1767–1772. [14] M. Ohtsu, K. Ono, A generalized theory of acoustic emission and Green’s functions in a half space, Journal of Acoustic Emission 3 (1984) 124–133. [15] A. Carpinteri, G. Lacidogna (Eds.), Earthquakes and Acoustic Emission, Taylor & Francis, Inc., ISBN 9780203936115, 2007. ¨ [16] W. Hufenbach, H. Richter, A. Langkamp, R. Bohm, Application of acoustic emission analysis for damage investigations in fibre and textile reinforced composites, in: Proceedings of the Conference on Damage in Composite Materials: Non Destructive Testing and Simulation (CDCM06), Stuttgart, 18.9.–19.9, 2006. [17] M. Ohtsu, Acoustic emission theory for moment tensor analysis, Research in Nondestructive Evaluation 6 (3) (1995) 169–184. [18] T. Shiotani, M. Ohtsu, K. Ikeda, Detection and evaluation of AE waves due to rock deformation, Construction Building Material 15 (2001) 235–246. [19] Monitoring Method for Active Cracks in Concrete by Acoustic Emission, Japan Construction Material Standards JCMS-IIIB5706, Japan: The Federation of Construction Material Industries, 2003. [20] A. Anastassopoulos, T. Philippidis, Clustering methodology for evaluation of acoustic emission from composites, Journal of Acoustic Emission 13 (1994) 11–22.