Mechanical Systems and Signal Processing 60-61 (2015) 971–985
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Determination of the probability zone for acoustic emission source location in cylindrical shell structures Ehsan Dehghan Niri, Alireza Farhidzadeh, Salvatore Salamone n Smart Structures Research Laboratory, University at Buffalo, The State University of New York, Buffalo, NY, USA
a r t i c l e i n f o
abstract
Article history: Received 24 October 2013 Received in revised form 1 July 2014 Accepted 9 February 2015 Available online 6 March 2015
This paper presents a probabilistic framework for acoustic emission (AE) source localization in cylindrical structures. Specifically, an approach based on unscented transformation (UT) is proposed to take into account uncertainty in time of flight measurements and wave velocity and eventually estimate AE source locations together with quantitative measures of confidence associated with those estimates. Experiments are carried out on a steel pipe instrumented with an array of six embedded piezoelectric disks. Results are compared with Monte Carlo simulations by the Kullback–Leibler divergence. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Acoustic emission Source localization Cylindrical shell structures Unscented transformation
1. Introduction Engineering structures such as pipelines, pressure vessels, and aircrafts are made of thin-walled shells whose performance and functionality is essential. Presence of defects, such as fatigue cracks, corrosion, impacts, in shell structures may severely compromise the overall soundness of these structures. Therefore, their proper assessment is crucial. To minimize the maintenance costs and to increase the operation lifetime, researchers and practitioners are increasingly interested in building advanced structural health monitoring (SHM) strategies [1]. In particular, SHM strategies based on acoustic emissions (AE) have been widely used [2–9] in the last three decades. In general, AEs are stress waves produced by sudden strain releases due to internal fractures [10]. A burst of energy can be released in the form of high-frequency sound waves from propagating cracks or from plastic deformations if a material is overstressed. Conventionally, an array of piezoelectric transducers is attached to the structure of interest to detect these waves and triangulation techniques are used to perform the critical task of damage localization [11–15]. These techniques work very well when the wave velocity (V) in the test material and the arrival time (ti) of the wave at all sensor locations is known. For example, considering an array of three sensors on the surface of a thin cylinder, the AE source can be identified by drawing three circles of radii (Ri), whose centers coincide with the three sensor locations, as shown in Fig. 1a. The radius Ri is the shortest path connecting the AE source and each sensor and it can be obtained by multiplying the time of arrival of the wave (ti) with the wave velocity (V). However, in real life, time of arrival measurements and wave velocity are in general affected by random and systematic errors. Random errors are caused by unknown and unpredictable changes in measurements, e.g. instrumentation noise. Systematic errors are mostly caused by the digital signal processing technique used for analyzing the time waveforms [16–18]. Therefore, rather than the damage being located at
n
Correspondence to: State University of New York at Buffalo, 212 Ketter Hall, Buffalo, NY 14260, USA. Tel.: þ 1 716 645 1523. E-mail address:
[email protected] (S. Salamone).
http://dx.doi.org/10.1016/j.ymssp.2015.02.004 0888-3270/& 2015 Elsevier Ltd. All rights reserved.
972
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
S1
AE source
S3 S2
S1
S2
S3
Fig. 1. AE Source localization using triangulation: (a) without uncertainty, (b) with uncertainty.
a single point (i.e., at the intersection of the circles as shown in Fig. 1a), the damage can be located anywhere in the dark overlapped region, as shown in Fig. 1b or it can remain undetected until failure of the structural system [19–21]. The ring’s width represents the uncertainty in the measured distance as a result of measurements uncertainty (i.e., ti and V). To consider these uncertainties, a probabilistic framework based on nonlinear Kalman Filtering techniques has been recently proposed by the authors for AE source localization in plate-like structures [16,19,20]. This framework is here extended to cylindrical shell structures. Specifically, an approach based on unscented transformation (UT) is proposed to take into account uncertainty in time of flight measurements and wave velocity and eventually estimate AE source locations together with quantitative measures of confidence associated with those estimates. 2. AE source localization in cylindrical shell structures Let’s consider a cylindrical structure instrumented with a sensor, and an AE source located on its surface, as shown in Fig. 2. It is known, that circumferential waves can propagate from the same source to the same sensor through infinite helical paths, as shown in Fig. 2a [13,22]. For structures with small wall thickness-to-diameter-ratio, these circumferential waves can be considered similar to Lamb waves in plates, by replacing the cylindrical structure with an equivalent “unwrapped” two dimensional plate [23]. Using the “unwrapped” representation, helical paths can be considered as a single AE source detected by multiple “virtual” sensors [13]. Assuming an arbitrary Cartesian coordinate system, the AE source is at the unknown coordinates (xs, ys), the sensor is located at (xi, yi), and the virtual sensors are placed at the vertically repeating positions, as shown in Fig. 2b. The distance between each pair of virtual sensors is πD, i 0 i where D is the diameter of the cylinder [13]. Also indicated in Fig. 2, three possible helical paths, whose lengths, l0 ; l1 i; l1 , can be defined as: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i ð1aÞ l0 ¼ ðxs xi Þ2 þ ðys yi Þ2 0
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxs xi Þ2 þ ðys ðyi þπDÞÞ2
ð1bÞ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxs xi Þ2 þ ðys ðyi πDÞÞ2
ð1cÞ
l1 i ¼ i
l1 ¼
Conventionally, AE source localization is performed using time-of-flight (t) measurements taken at multiple receiving points. Particularly, given an array of n sensors, if t m is the travel time of the first arrival waveform to reach the first triggered sensor (master sensor), and Δt mi is the time difference between master sensor and the ith sensor, the following equations can be obtained: m
ð2Þ
i
ð3Þ
lmin t m V ¼ 0 lmin ðt m þΔt mi ÞV ¼ 0
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
AE Source
973
Theoretical Cut Line
Theoretical Cut Line
Virtual sensor
AE Source
Theoretical Cut Line
Virtual sensor
Fig. 2. The first three shortest helical paths (a) 3-D view, (b) unwrapped cylinder (Cartesian coordinate). i
m
where V is the wave velocity and lmin and lmin are the shortest paths (or Geodesic) between AE source and, ith sensor and master sensor, respectively. In a plane structure, the shortest path is a straight line connecting the AE source and the sensor [11]. However, on a curved structure the shortest path should be determined by writing the equation for the length of a curve and then minimizing this length using the calculus of variations [24,25]. For example, for the case shown in Fig. 2, the shortest path between the AE source and the ith sensor should be defined as: i
0i
i
i
lmin ¼ minðl0 ; l1 ; l1 Þ
ð4Þ
Combining Eq. (2) and Eq. (3) the following set of (n 1) nonlinear equations can be obtained: i
Δt mi ¼
m
lmin lmin V
ð5Þ
It is well known that given the wave velocity, at least three sensors are necessary to locate an AE source in plate-like structures. However, this may lead to ambiguous answers for cylindrical structures [14]. The addition of a forth sensor located unsymmetrically with respect to the longitudinal axis could be used to overcome this issue [14]. The additional sensor results in a system of nonlinear equations which is over-determined, that is, the number of equations, (i.e., n 1), are more than the number of unknowns, (i.e., xs, ys). Furthermore, since the time difference measurements are in general affected by noise, Eq. (5) is not satisfied exactly. Therefore, the problem of finding the AE source location coordinates (xs, ys) could be phrased as an optimization problem where the objective is to minimize the weighted square of residual errors while simultaneously satisfying the constraint defined by Eq. (4). This can be expressed mathematically as: T minimize J xs ; ys ¼ Z~ Z W Z~ Z i i 0i i subject tolmin ¼ min l0 ; l1 ; l1
ð6Þ
h i f mi where Z~ is a vector which contains the noisy time difference measurements, (i.e. Z~ ¼ Δt
ðn 1Þ1
), Z is a vector which
components are defined in Eq. (5), (i.e. Z ¼ ½Δt mi ðn 1Þ1 ), and W is a weighting matrix defined as the inverse of the error covariance matrix: 2
σ2 e m1 6 Δt 6 W ¼ R1 ¼ 6 ⋮ 4 0
⋯
0
⋱
⋮
⋯
σ2 e mðn 1Þ Δt
31 7 7 7 5 ðn 1Þðn 1Þ
ð7Þ
Sensor 2
Signal 2
Sensor n
Signal n
y
AE Source Location Estimation
Sensor 1
Evaluating the Objective Function With Uncertainty Considerations J
Signal 1
Theoretical Cut Line
Mesh Gridding and Evaluating Time of Flight at each grid point
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
Time of Flight Measurements
974
x Unwrapped Cartesian Theoretical Cut Line Coordinates
Fig. 3. AE source localization algorithm flowchart.
If we assume that the arrival time (ti) are mutually independent Gaussian random variables with variance σ 2ti the variance f mi can be defined as: of the time difference Δt σ 2g ¼ σ 2ti þσ 2tm Δt mi
ð8Þ
The optimization problem defined in Eq. (6) is in general non-convex. A physical intuition behind this problem is due to the fact that because of the geometry of cylinder in some cases circles shown in Fig. 1 might have two intersections (a circle can also intersect with itself). Several algorithms have been developed in literature for solving this kind of problems, such as Genetic Algorithms (GA) [26,27], Particle Swarm Optimization (PSO) [28], and Mesh Grid Optimization (MGO) [29]. In this paper,a MGO algorithm was used, with 1 mm resolution in both x and y directions. Therefore, the AE source coordinates xs ; ys was estimated by evaluating the objective function in Eq. (6) at gridded locations on the unwrapped surface. The resulting absolute minimum value represented the unknown AE source location. The overall source localization algorithm proposed here is shown in Fig. 3. 3. Time of flight measurement and its systematic error Several methods can be used for time of flight measurement of acoustic signals, including, threshold crossing, crosscorrelation, wavelet transform, and curve fitting approaches [10]. In this work, threshold crossing of the signal envelope was used to calculate the time of flight of each waveform. The signal envelope can be calculated using the analytic representation of a real-valued signal that is: xA ðt Þ ¼ xðt Þ þ ixH ðt Þ where xðt Þ is a real valued signal, i is the imaginary unit, and xH ðt Þ is the Hilbert transform of xðt Þ defined as: Z 1 þ1 1 dτ xðτÞ x H ðt Þ ¼ π 1 t τ The envelope can be defined as: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Aðt Þ ¼ xðt Þ2 þ xH ðt Þ2
ð9Þ
ð10Þ
ð11Þ
The time instant in which the envelope of the signal recorded by the ith sensor, reaches the threshold (At), is considered as the time of flight or arrival time t i . The sampling frequency ðf s Þ used in the data acquisition system is in general considered as a source of uncertainty [30]. This is a systematic uncertainty that can be determined analytically as a function of the sampling frequency [10,16,31]. Let’s consider a discretized signal xðndt Þ, where dt is the sampling period, dt ¼ 1=f s , and n is an integer number. Similarly the analytic signal Aðt Þ can be discretized as Aðndt Þ. In general, the exact time of flight ðt i Þ of the ith sensor can be expressed as: t i ¼ ðk þ λÞdt
ð12Þ
where k is an integer and λ is a real number in the interval [1 0], while the time of flight t~ i measured by the threshold crossing method is: t~ i ¼ kdt
ð13Þ
Therefore a systematic error εt i can be defined as: εti ¼ t i t~ i ¼ λdt For illustrative purpose, Fig. 4 shows the error for a discretized signal.
ð14Þ
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
975
Amplitude [Volt]
Threshold
Time [msec] Fig. 4. Sampling error due to discretization.
The time of arrival variance is directly related to this systematic error λ so that its variance σ 2t should be identified according to this error. In this paper we assumed time of flight systematic error as zero mean Gaussian noise with standard deviation σ t equal to the uniform distribution band, i.e. σ t ¼ dt ¼ ð1=f s Þ. Under this assumption the variance of time of flight differences error σ 2 in Eq. (7) can be defined as: e mi Δt 2 1 σ 2e ¼ σ 2ti þ σ 2tm ¼ 2 ð15Þ Δt mi fs It should be noted that, in this paper just the systematic error due to the digital time resolution was considered in the time of flight measurement. However, if other sources of uncertainty such as random noise, were known, their effect on the variance of time of flight differences error σ 2 should be taken into account. e mi Δt 4. Uncertainty propagation using unscented transformation Let us consider a nonlinear function f( ) that relates a random variable Y to an N-dimensional random variable X as: Y ¼ f ðX Þ
ð16Þ
Assume X has mean X and covariance PX. We want to estimate mean and covariance of Y. The simplest method to transform uncertain parameters and estimate their mean and covariance is the Monte Carlo simulation. The Monte Carlo simulation random sampling. Particularly, to capture mean Y and covariance PY a large number of relies on repeated
samples X 1 ; X 2 ; …X Nm , are randomly drawn from X and the function f( ) is evaluated for each sample. Transferring these
samples through the function f( ), we will have a set of N m transferred samples Y 1 ; Y 2 ; Y 3 ; …; Y Nm . The mean Y and covariance P Y are estimated as: Y¼
Nm 1 X Yi Nm 1
PY ¼
ð17aÞ
Nm T 1 X Y i Y Y i Y Nm 1
ð17bÞ
The main drawback of this method is that it requires a large number of calculations and can be prohibitive when each calculation involves a long computational cost. In this paper an alternative approach based on unscented transformation (UT) was used. The UT is a method for calculating the statistics of a random variable which undergoes a nonlinear transformation [32,33]. If the number of random variable included in the input vector X is N, to calculate the statistics of Y (mean Y and covariance P Y ) a matrix of 2N þ1 samples called sigma points χ i and their associate weights wi are “deterministically” chosen so that the first two moments of these points match those of the prior distribution of X. The sigma points and their corresponding weights are defined as follows [32,33]: χ ð0Þ ¼ X;
κ w0 ¼ N þ κ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðiÞ ðN þκ ÞP x ; χ ¼ Xþ
wi ¼
i
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi χ ði þ N Þ ¼ X ðN þκ ÞP x ; i
1 ; 2ðN þ κÞ
wi þ N ¼
i ¼ 1…N
1 ; 2ðN þκ Þ
i ¼ 1…N
ð18Þ
976
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
where N isthe number of uncertain parameters or size of X, κ is a real number that determines the spread of sigma points pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðN þ κÞP x is the ith row or column of the matrix square root of ðN þ κ ÞPx and wi is the weight associate with around X, i the ith sigma point. κ is any positive or negative number providing N þ κ a0; κ ¼ 0 was used in this paper. The transformation of each sigma point through the nonlinear function f( ), Eq. (16), gives: ð19Þ Ψi ¼ f χ i The mean and covariance for Y can be approximated by the weighted average of transformed points as: Y
2N X
wi Ψi
ð20Þ
0
Py
2N X
T wi ðΨi YÞ Ψi Y
ð21Þ
0
The UT has two main advantages over Monte Carlo simulation: (1) only a small number of sigma points deterministically chosen, is required therefore, no random sampling is needed; whereas; (2) it is faster than Monte Carlo. 5. Probability zone determination of AE source location This section presents a procedure to determine a probability zone for AE source location in cylindrical structures. First, the probability distribution of the AE location coordinates xs ; ys as a result of uncertain parameters is estimated. Then, we predict shape and location of a zone within which point source has a certain probability of existence. 5.1. AE source location distribution Traditionally, AE source localization is performed assuming the source location, to be a deterministic quantity. Instead of solving for the point estimate for source location, we are interested in probability distribution for their value due to uncertainty in time of flight and wave velocity. Thus, the location of damage is assumed to be a random vector with total uncertainty characterized by the probability distribution function (pdf). If we assume that the source location xs ; ys have a bivariate Gaussian distribution, only its mean xs ; ys and covariance matrix P are needed. It should be noted that the source localization algorithm described in Section 2 can be interpreted asa nonlinear transformation process that resembles the nonlinear function f( ) in Eq. (16), in which Y ¼ xs ; ys , Y ¼ xs ; ys and P y ¼ P. This algorithm processes some uncertain parameters, say input X, and estimate the source location coordinates xs ; ys . The input X can include both or either of time f mi and wave velocity V. If uncertainty in both time of flight and wave velocity is considered the input of flight differences Δt X and its covariance matrix P X are defined as: i 2h 3 f mi Δt ðn 1Þ1 5 ð22aÞ X ¼4 V n1 " PX ¼
R
0
0
σ 2V
# ð22bÞ
where R is defined in Eq. (7), V and σ 2V are mean and variance of the wave velocity distribution. Given the input X and its covariance P X sampling based approaches, such as Monte Carlo and UT, can be used to find the pdf of the AE source location. The number of parameters N included in X is n that is (n 1) time of flight differences and the wave velocity. Fig. 5 schematically shows the propagation of uncertain parameters through the AE source localization algorithm, using Monte Carlo and UT, respectively.
In the Monte Carlo simulation, as it was discussed, a set of Nm sample points, X 1 ; X 2 ; …; X Nm , was randomly selected and a set locations, fðxs ; ys Þ1 ; ðxs ; ys Þ2 ; …; ðxs ; ys ÞNm g, was estimated using the algorithm described in Section 2. Then, of source mean xs ; ys and covariance P were calculated using Eqs. (17a) and (17b): " # " # Nm xs xs 1 X ¼ ð23aÞ ys y Nm 1 s MC
P MC ¼
Nm 1 X Nm 1
i
"
xs ys
#
" i
xs ys
#! "
xs ys
#
" i
xs ys
#!T ð23bÞ
The subscript “MC” was used to indicate that these values were obtained from a Monte Carlo simulation. Conventionally a sensitivity test for determining the minimum number of required Monte Carlo samples N m , is carried out. This test increases the number of samples till the stability in mean and covariance is achieved. The number of samples in this paper was set to 1000.
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
977
Fig. 5. Uncertainty transformation using sampling approach: (a) Monte Carlo, and (b) unscented transformation.
The UT approach, as it said earlier, addresses the uncertainty propagation of X by using a deterministic and efficient
sampling approach. Specifically a set of ð2n þ1Þ points called sigma points, χ ¼ χ ð0Þ ; χ ð1Þ ; χ ð2Þ ; …; χ ð2nÞ , selected using (" # " # xs xs ; ; Eq. (18), were propagated using the source localization algorithm, to generate transformed points, Ψ ¼ ys ys 0 1 " # " # xs xs ; …; . The mean and covariance of the estimated location were calculated using Eqs. (20) and (21): ys ys 2
2n
"
xs ys
#
wi Ψi
ð24Þ
0
UT
P UT
2n X
2n X
" wi Ψi
xs
#!
"
xs
#!T
ð25Þ Ψi ys ys The subscript “UT” was used to indicate that these values were obtained by an unscented transformation approach. Having the mean and covariance of estimated AE source location, the next step is to find a region in which with the certain level of confidence the AE source is located. 0
5.2. Determination of the covariance ellipse with α% confidence interval In this section, using the mean and covariance of the estimated AE source location distribution, we wish to determine a region such that the probability of the actual point laying within this region is α. In general a bivariate AE source Gaussian distribution can be expressed as: " # " # !! xs xs 1 p xs ; ys ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiexp Δ ; ;P ð26Þ ys ys 2 ð2π Þ jP j
978
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
" where " from
#
xs
" and P are mean and covariance of AE source distribution. The Δ
ys #
xs
" to mean
ys " Δ
# "
xs
;
ys
xs
xs
#
ys
ys
"" ¼
xs
#
"
ys
xs
##T
ys
P 1
""
xs ys
#
"
xs
γ
xs ys
#
"
xs
##T
ys
;
xs
#
ys
;P
! is the Mahalanobis distance
##
ys
A particular Mahalanobis distance γ, gives an ellipsoid centered at the mean ""
# "
under covariance P, defined as:
ys # ! ;P
xs
P
1
""
xs ys
#
"
xs ys
ð27Þ "
xs ys
# as:
## ¼0
ð28Þ
The confidence interval α associated with this ellipse is the probability mass enclosed by it. In the other word, a particular confidence interval α uniquely identifies an ellipse with a fixed Mahalanobis distance γ [34]. Thus we wish to find γ associated with α confidence interval to identify the region for AE source as a probability zone or covariance ellipse in the form of Eq. (28) with α% confidence interval. Since the Mahalanobis distance in Eq. (28) is sum of two Gaussian random variables, the Mahalanobis distance is distributed according to the chi-squared distribution with two degrees of freedom. Therefore, the Mahalanobis distance γ associated with α% confidence interval can be determined by the inverse cumulative distribution function of two degree chi-squared distribution up to confidence interval α as: α¼
Z
γ 0
e2 dt 2 t
ð29Þ
Solving γ for given α we have: γ ¼ 2 lnð1 αÞ
ð30Þ
Having γ the equation for the covariance ellipse with the α% confidence interval can be drawn from Eq. (28). More references for defining covariance ellipse can be found in [35].
5.3. Kullback–Leibler divergence for quantitative evaluation of UT method In order to quantitatively evaluate the AE distribution obtained by UT and Monte Carlo, the Kullback–Leibler divergence, (DKL ) was used. In probability theory the DKL is a measure of the difference between two probability distributions. The DKL for two continues probability distributions pðxÞ and qðxÞ is defined as [36]: DKL ðp J qÞ ¼
Z
þ1
pðxÞln 1
pðxÞ dx qðxÞ
ð31Þ
Circular PZT
Fig. 6. PZT disk.
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
979
where x is a random variable. The DKL is a non-symmetric and non-negative value. DKL is zero if and only if pðxÞ ¼ qðxÞ in almost everywhere. The DKL between two bivariate Gaussian distributions resulted from Monte Carlo and UT can be calculated as: " # " # ! " # " # !T " # " # ! xs xs xs xs xs xs 1 1 ¼ 12ðtr P MC J P UT þ P MC DKL ys ys ys ys ys ys ð32Þ UT
MC
detðP UT Þ 2 lnðdetðP ÞÞ " # " MC Þ # xs xs , where ys ys UT
MC
UT
MC
UT
,P UT , and P MC are defined in Eqs. (23a, 23b) and (24).
MC
6. Experimental setup Experiments were carried out using a 3000 mm long steel pipe with 152.4 mm diameter and 2.1 mm thickness, to validate the proposed algorithm. The pipe was instrumented with an array of six piezoelectric disks (see Fig. 6) permanently attached to the surface of the pipe using a Loctite instant adhesive. The specimen was gridded longitudinally and
Group 1 Group 3 S1 Group 2
S6
S3 60o
S5
S4 S2
Group 1 (S1, S2)
S4 S2
1m Group 2 (S3, S4)
S1 S3
Amplifiers
Data Acquisition
S6 S5
1m
D
Group 3 (S5, S6)
Fig. 7. Experimental setup: (a) sensor layout, (b) test specimen and data acquisition system.
Fig. 8. Wave velocity distribution of experimental data and fitted Gaussian density function.
980
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
S1
Theoretical Cut Line
1
4
z [cm]
3
2
S3 6
7
10
S6
9
5
8
S5
11
S2
Fig. 9. Simulated AE sources: (a) 3D view, (b) 2D view.
Table 1 Sensor locations. i
1
2
3
4
5
6
xi [cm] yi [cm]
200 5
200 28.94
100 12.98
100 36.92
0 20.96
0 44.9
circumferentially with 40 mm and 1/8 of pipe perimeter increment respectively, as shown in Fig. 7. Each sensor was connected to a 40 dB preamplifier. The AE system included an eight-channel high-speed data acquisition board (Physical Acoustics Corporation Micro-II PAC) with sampling frequency fs ¼3 MHz and dedicated software (AEwin) for signal processing and storage. A threshold was set 6 dB above the background noise for each sensor. Post-processing of the AE signals was performed with a PC running a Matlab software code implemented by the authors. The wave velocity was estimated by performing pencil lead breaks at known locations, and measuring the time of flight differences Δt ij between two sensors, that is: V¼
li lj Δt ij
ð33Þ
where li and lj were the minimum distance between the simulated AE source and the ith and jth sensor, respectively. A total of 183 wave velocities were measured at different locations using different pair of sensors. The results are summarized in Fig. 8. Mean and standard deviation of the measured wave velocity were V ¼ 5277 m=s and σ V ¼ 244 m=s, respectively. The sensor layout is shown in Fig. 9. Table 1 lists the sensors’ coordinates in the unwrapped coordinate system. 7. Results To validate the proposed approach, AE sources were generated by pencil-lead breaks at eleven grid locations (see Fig. 9). At each location, the first four triggered sensors were used as input to the proposed algorithm. A typical objective function J [ ] is plotted in Fig. 10. Specifically, this figure shows the objective function defined in Eq. (6) for an AE source generated at the location 2, (see Fig. 9). Clearly, this function is non-convex with a local minimum. Therefore, conventional gradientbased optimization methods, which tend to get stuck in local minima, are not suitable for TOF-based source localization methods in cylindrical structures. The data used to populate the input vector X and its covariance matrix P X , including time of flights measurements and wave velocity are listed Table 2. Fig. 11 shows the probability zone/covariance ellipse corresponding to α ¼ 90% confidence interval, estimated using a Monte Carlo simulation. A set of 1000 sample points fX 1 ; X 2 ; …; X 1000 g was randomly selected
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
981
Fig. 10. Weighted square of residual error for simulated point 2 and 1 mm resolution in x and y direction.
Table 2 Extracted time of arrivals for simulated point 2 and the first four triggered sensors. Data
m=s
t 2 ½ms
t 1 ½ms
t 3 ½ms
t 4 ½ms
V
0.2557
0.2693
0.3813
0.3827
5277
σ V ½m=s
f s ½MHz
244
3
Fig. 11. Estimated locations of point 2 for 1000 Monte Carlo samples.
from the input vector X and its covariance matrix P X , and a set of source locations xs ; ys 1 ; xs ; ys 2 ; …; xs ; ys 1000 , was estimated using the localization algorithm described in Section 2. From this data set, mean and covariance were calculated using Eqs. (23a) and (23b). These values are named “true” since the number of samples was large enough to accurately estimate the true value of mean and covariance. The axes are the two AE source coordinates of interest. In this figure, the mean-value and the covariance ellipsoid computed by the Monte Carlo simulations are indicated. The individual realizations of the Monte-Carlo simulations are also plotted to show the actual dispersion of the estimated points. The results obtained using the UT are shown in Fig. 12. Since the input vector X contains four uncertain parameters, the minimum number of sigma points required is nine (i.e.,2 4 þ 1). Table 3 lists the selected sigma points along with the corresponding weights determined using Eq. (18). Then, each sigma point was propagated using the source localization algorithm described in Section 2, to generate the transformed points listed in Table 4. Finally, mean and covariance of the estimated AE source location were calculated using Eq. (24) and Eq. (25). In Fig. 12, the mean-value and the covariance ellipsoid computed by the UT approach are indicated. Also plotted are the transferred sigma points.
982
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
Fig. 12. Estimated locations of point 2 for 9 sigma points samples using UT. Table 3 Generated sigma points for simulated point 2. Sigma point
χ ð0Þ
χ ð1Þ
χ ð 2Þ
χ ð3Þ
χ ð4Þ
χ ð5Þ
χ ð 6Þ
χ ð7Þ
χ ð8Þ
f 21 [ms] Δt f 23 [ms] Δt f 24 [ms] Δt
0.01366
0.014603
0.01366
0.01366
0.01366
0.012717
0.01366
0.01366
0.01366
0.12566
0.12566
0.12660
0.12566
0.12566
0.12566
0.12472
0.12566
0.12566
0.1270
0.1270
0.1270
0.12794
0.1270
0.1270
0.1270
0.12606
0.1270
V [m/s] wi
5277
5277
5277
5277
5763
5277
5277
5277
4791
0
1=2 4
1=2 4
1=2 4
1=2 4
1=2 4
1=2 4
1=2 4
1=2 4
Table 4 Transferred sigma points for simulated point 2.
xs [cm] ys [cm]
Ψ0
Ψ1
Ψ2
Ψ3
Ψ4
Ψ5
Ψ6
Ψ7
Ψ8
183.39 23.24
183.29 23.74
183.49 23.34
183.49 23.24
186.59 23.04
183.39 22.94
183.29 23.24
183.19 23.44
180.19 23.44
Fig. 13. Estimated and true (Monte Carlo) mean and covariance ellipse for 90% confidence interval.
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
983
A comparison between the probability zone estimated using the UT and Monte Carlo simulations, is shown in Fig. 13. It can be seen that the mean value and covariance computed by the UT are close to the “true” value computed by the Monte Carlo simulation.
S1
Theoretical Cut Line
1
4
S3
7 10
S6 z [cm]
3
6
9
2 5
S2
S1 4
S3
10
S6
6
3
z [cm]
Theoretical Cut Line
7
1
S5
11
8
9
2 5
8
11
S5
S2
Fig. 14. Probability zone for all 12 simulated points (a) Monte Carlo, (b) UT.
Fig. 15. Probability zone estimated for all 12 simulated points in the unwrapped Cartesian coordinate system (a) Monte Carlo, (b) UT.
984
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
Fig. 16. Kullback–Leibler divergence of AE source distributions of UT and Monte Carlo for 11 simulated points.
The true and estimated probability zones with 90% confidence interval for the AE sources generated at the 11 locations are depicted in Fig. 14. It can be observed that all the probability zones encompass the actual source location and the approximated mean calculated using UT are very close to the actual locations. For the sake of clarity the results in the unwrapped cylinder representation are shown Fig. 15. In order to quantitatively evaluate the performance of UT and Monte Carlo, the DKL is calculated using Eq. (32). The results are shown in Fig. 16 for the 11 simulated AE sources. The small values of DKL indicates that UT can estimate with a good accuracy the AE location distribution. To have better idea about the DKL values presented in Fig. 16, it is worthy to compare the DKL result of point 2 to other points. The DKL for point 2 is 0.018 and the Monte Carlo and UT probability zones corresponding to this value are compared in Fig. 13.
8. Conclusions Traditionally, AE source localization is performed assuming the source location, to be a deterministic quantity. Instead of solving for the point estimate for source location, we were interested in probability distribution for their value due to uncertainty in time of flight and wave velocity. Thus, the source location was assumed to be a random vector with total uncertainty characterized by the probability distribution function. In this paper an approach based on unscented transformation (UT) was presented to estimate the mean and covariance of the AE source location distribution. To evaluate the performance of UT, Monte Carlo simulations were carried out, and a Kullback–Leibler divergence was used. It was shown that UT was able to accurately estimate the mean and covariance of the AE location using a set of samples smaller than the Monte Carlo simulations. The accuracy the proposed framework demonstrated its potential application in real-time SHM applications. However, more formal tests need to be conducted to verify the robustness of the approach to more practical defect and damage types, such as cracks, corrosion, etc.
Acknowledgments Funding provided by the NYS Pollution Prevention Institute through a grant from the NYS Department of Environmental Conservation. Any opinions, findings, conclusions or recommendations expressed are those of the author(s) and do not necessarily reflect the views of the Department of Environmental Conservation. References [1] F.J. Weisweiler, G.N. Sergeev, Non-destructive Testing of Large-diameter Pipe for Oil and Gas Transmission Lines, VCH, Wienheim, 1987. [2] H.S. Jee, J.O. Lee, N.H. Ju, C.H. So, J.K. Lee, Acoustic emission source location of type-II composite pressure vessel with artificial defect. Sixth NDT in Progress 2011, International Workshop of NDT Experts, Prague, 2011. [3] H.S. Jee, J.O. Lee, Acoustic emission of composite vessel, in: N. Hu (Ed.), Composites and Their Applications, InTech, 2012, pp. 61–90. [4] P. Nivesrangsan, J.A. Steel, R.L. Reuben, Source location of acoustic emission in diesel engines, Mech. Syst. Sig. Proc. 21 (2007) 1103–1114. [5] H. Chen, Y. Lu, H. Zhao, Y. Sun, P. Wang, J. Huang, Automatic gauge control hydraulic cylinder state identification using modified image based acoustic emission profile, in: Proceedings of the Second International Conference on Computer Science and Electronics Engineering, 2013, pp.2153–2159. [6] C. Jomdecha, A. Prateepasen, P. Kaewtrakulpong, Study on source location using an acoustic emission system for various corrosion types, NDT E Int. 40 (2007) 584–593.
E. Dehghan Niri et al. / Mechanical Systems and Signal Processing 60-61 (2015) 971–985
985
[7] L. Gaul, S. Hurlebaus, L. Jacobs, Localization of a synthetic acoustic emission source on the surface of a fatigue specimen, Res. Nondestructive Eval. 13 (2001) 105–117. [8] T. Kundu, S. Das, K.V. Jata, Detection of the point of impact on a stiffened plate by the acoustic emission technique, Smart Mater. Struct. 18 (2009) 035006. [9] R.K. Miller, A.A. Pollock, D.J. Watts, J.M. Carlyle, A.N. Tafuri Jr, JJY., A reference standard for the development of acoustic emission pipeline leak detection techniques, NDT E Int. 32 (1999) 1–8. [10] J.A. Baron, S.P. Ying, Acoustic emission source location, in: K. Miller, P. Mclntire (Eds.), Nondestructive Testing Handbook, American Society for Nondestructive Testing, 1987, pp. 135–155. [11] Tobias. Acoustic emission source location in two dimensions by an array of three sensors, Int. J. Nondestr. Test. 9 (1976) 9–12. [12] M. Asty, Acoustic emission source location on a spherical or plane surface, NDT E Int. 11 (1978) 223–226. [13] D.J. Yoon, YH. Kim, OY. Kwon, New algorithm for acoustic emission source location in cylindrical structures, J. Acoust. Emission 9 (1992) 237–242. [14] P. Barat, P. Kalyanasundaram, B. Raj, Acoustic emission source location on a cylindrical surface, NDT E Int. 26 (1993) 295–297. [15] E. Lympertos, E. Dermatas, Best sensors position for accurate location of acoustic emission sources on cylindrical surfaces, in: International Conference on Computational & Experimental Engineering and Science, Corfu Hellas, 2003. [16] E. Dehghan Niri, S. Salamone, A probabilistic framework for acoustic emission (AE) source localization in plate-like structures, Smart Mater. Struct. 21 (2012) 035009. [17] W.F. Hartman, J.W. McElroy, Acoustic Emission Monitoring of Pressurized Systems, American Scosiety for Testing Materials, 1979. [18] W.H. Sachse, S. Sancar, Acoustic Emission Source Location on Plate-Like structures using a Small Array of Transducers. United State Patent 4592034, 1986. [19] E. Dehghan Niri, A. Farhidzadeh, S. Salamone, Adaptive multisensor data fusion for acoustic emission source localization in noisy environment, Int. J. Struct. Health Monit. 12 (2013) 59–77. [20] E. Dehghan Niri, A. Farhidzadeh, S. Salamone, Nonlinear Kalman filtering for acoustic emission source localization in anisotropic panels, Ultrasonics (2013), http://dx.doi.org/10.1016/j.ultras.2013.07.016. [21] A. Farhidzadeh, E. Dehghan-Niri, S. Salamone, B. Luna, A. Whittaker, Monitoring crack propagation in reinforced concrete shear walls by acoustic emission, J. Struct. Eng. ASCE (2012), http://dx.doi.org/10.1061/(ASCE)ST.1943-541X.0000781. [22] K.R. Leonard, M.K. Hinders, Guided wave helical ultrasonic tomography of pipes, J. Acoust. Soc. Am. 114 (2003) 767–774. [23] J. Li, J.L. Rose, Natural beam focusing of non-axisymmetric guided waves in large-diameter pipes, Ultrasonics 44 (2006) 35–45. [24] R. Gangadharan, G. Prasanna, M.R. Bhat, C.R.L. Murthy, S. Gopalakrishnan, Acoustic emission source location and damage detection in a metallic structure using a graph-theory-based geodesic approach, Smart Mater. Struct. 18 (2009) 115022. [25] R. Gangadharan, G. Prasanna, M.R. Bhat, C.R.L. Murthy, S. Gopalakrishnan, Acoustic emission source location in composite structure by Voronoi construction using geodesic curve evolution, J. Acoust. Soc. Am. 126 (2009) 2324–2330. [26] E. Dehghan Niri, S.M. Zahrai, A. Mohtat, Effectiveness-robustness objectives in MTMD system design : an evolutionary optimal design methodology, Struct. Control Health Monit. 17 (2010) 218–236. [27] P.T. Coverley, W. Staszewski, Impact damage location in composite structures using optimized sensor triangulation procedure, Smart Mater. Struct. 12 (2003) 795–803. [28] I. Muller, F.K. Chang, D. Building, S. Mall, Model-based impact monitoring by inverse methods using particle swarm optimization, in: Proceedings of IMAC XXVII, 2009. [29] T. Kundu, S. Das, K.V. Jata, Health monitoring of a thermal protection system using Lamb waves, Int. J. Struct. Health Monit. 8 (2008) 29–45. [30] G. Andria, F. Attivissimo, N. Giaquinto, Digital signal processing techniques for accurate ultrasonic sensor measurement, Measurement 30 (2001) 105–114. [31] D. Marioli, C. Narduzzi, C. Offelli, D. Petri, E. Sardini, A. Taroni, Digital time-of-flight measurement for ultrasonic sensors, IEEE Trans. Instrum. Meas. 41 (1992) 93–97. [32] S.J. Julier, J.K. Uhlmann, New extension of the Kalman filter to nonlinear systems, in: Proceedings of AeroSens: 11th International Symposium on Aerospace/Defense Sensing Simulation and Controls, 1997, pp. 182–193. [33] S.J. Julier, J.K. Uhlmann, H.F. Durrant-whyte, A new method for the nonlinear transformation of means and covariances in filters and estimators, IEEE Trans. Autom. Control 45 (2000) 477–482. [34] V. Chew, Confidence, prediction, and tolerance regions for the multivariate normal distribution, J. Am. Stat. Assoc. 61 (1966) 605–617. [35] A. Alexandersson, Graphing confidence ellipses: an update of ellip for Stata 8, Stata J. (2004) 242–256. [36] S. Kullback, Information Theory and Statistics, Wiley Publications in Statistics, 1959.