Table of Contents Author Guidelines Submit a Manuscript
Shock and Vibration
Volume 2015, Article ID 964805, 10 pages
http://dx.doi.org/10.1155/2015/964805
Research Article

Performance Improvement of Ensemble Empirical Mode Decomposition for Roller Bearings Damage Detection

Dynamics & Identification Research Group, Department of Mechanical and Aerospace Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy

Received 10 October 2014; Revised 11 February 2015; Accepted 24 February 2015

Academic Editor: Ahmet S. Yigit

Copyright © 2015 Ali Akbar Tabrizi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Ensemble empirical mode decomposition (EEMD) is a noise assisted method widely used for roller bearing damage detection. However, to successfully handle this technique still remains a great challenge: identification of two effective parameters (the amplitude of added noise and the number of ensemble trials), which affect the performances of the EEMD. Although a number of algorithms or values have been proposed, there is no robust guide to select optimal amplitude and the ensemble trial number yet, especially for early damage detection. In this study, a reliable method is proposed to determine the suitable amplitude and the proper number of trials is investigated as well. It is shown that the proposed method (performance improved EEMD) achieves higher damage detection success rate and creates larger Margin than the original algorithm. It leads to a substantially low trial numbers required to achieve perfect labelling of samples; in turn this fact leads to considerably less computational cost. The number of real vibration signals is analysed to verify effectiveness and robustness of the proposed method in discriminating and separating the faulty conditions.

1. Introduction

Modern rotating machines become more precise and automatic, fast and costly. As an obvious consequence, their lifetime is extended as much as possible, and this fact implies a strong request of increasing reliability and capacity of detecting faults at a very early stage. Through the processing of collected vibration signals and the extraction of significant information, it is possible to detect even small defects appearing on bearings during their life. Several signal processing techniques exist to decompose a signal and extract informative features. Randall and Antoni have broadly treated the background of a number of successful diagnostic methods [1]. The EMD is another recent adaptive data driven technique [2], to decompose a multicomponent signal into several elementary intrinsic mode functions (IMFs) and has been widely applied to fault diagnosis of rotating machines. However, there exist some drawbacks of the method such as the stopping criterion for sifting process, the mode mixing, and the border effect problem. The intermittency of the detected extrema, which belong to the different orthogonal components, is the main reason of mode mixing effect. EEMD is a noise assisted data analysis method and has been recently proposed to eliminate the mode mixing problem of the EMD technique [3]. Essentially, the EEMD repeatedly decomposes the original signal with added white noise into a series of IMFs, by applying the original EMD process. The means of the corresponding IMFs during the repetitive process is considered as the final EEMD decomposition result. Since white noise is added throughout the entire signal decomposition process, mode mixing is effectively eliminated. The EEMD has been already used to detect rotating machine faults such as defective bearings and gears in the past few years [4].

However, another challenge still exists: how to better identify the two effective parameters (the amplitude of added noise and the number of ensemble trials), which affect the performance of the EEMD. If the amplitude of the added noise is too small relative to the original signal, a considerable mode mixing improvement cannot be achieved. On the other hand, if the amplitude of the added noise is too high, it will create some redundant IMF components which lead to misinterpretation of the analysis result. In addition, although an infinite number of ensemble trials are required to completely cancel out the effect of the added white noise, too many trial numbers would increase the computational cost. Wu and Huang [3] suggested the value of 0.2 of standard deviation of the original signal as the amplitude of the added white noise and a few hundred for trial number of ensemble. It has been shown in various cases that such an amplitude is not appropriate. Zhang et al. [5] suggested using a band-limited white noise to decrease the computational cost. Analysing a simulated signal, it was concluded that appropriate range of SNR (signal-to-noise ratio based on signal power) is (50–60) dB. However, they used another range ((0.01–0.1)), which is outside of the suggested SNR. A nonstationary signal was constructed to mimics realistic vibration signals measured from rolling bearing and the appropriate range of SNR was considered (49–58) dB for the vibration signals. Applying the EEMD to the simulated signals, it was obtained that when the number of ensemble trials is 100, the corresponding correlation coefficient approaches 0.95. Using the modified EEMD method, the acceptable results were achieved approximately after 70 ensemble trials, instead of the 100 trials suggested for the original EEMD method. For real data (acceleration signals), it was shown that the percentage improvement of the computational efficiency (the consumed time ratio) varies from 30% to 45%, depending on the operating conditions. Guo and Tse [6] investigated the influence of the parameters setting on the results of reducing the mode mixing problem using a simulated signal. The effects of frequency and amplitude ratio of two different parts of the simulated signal (the high frequency and low frequency components) were investigated as well. The investigated amplitudes were considering again coefficients of standard deviation of the original signal (0.01, 0.1, 0.2, and 0.3). As real data is noisy (produced by other industrial equipment) and the amplitudes and composition of frequency are unknown, lower amplitude of noise was added and more number of ensemble trials applied (0.1 of standard deviation of the original signal for amplitude and 3000 for ensemble trial number). As only one specific operating condition with a single predefined amplitude was investigated, it would not represent a reliable guideline for properly setting the best parameters for real signals. Lin [7] tried to provide guidance on choosing the appropriate amplitude and reduce the tremendous time waste occurring in the EEMD method. An optimal interval was suggested that lies between the square root of the average power of the weak sinusoid component and that of the weak transient component. When the amplitude is selected from the mentioned interval, Pearson’s correlation coefficients of the components reach their maximum value. Taking into consideration that only one specific gearbox vibration signal was investigated to verify the suggested procedure, its performance does not seem too reliable to identify small defects. Furthermore, it seems difficult to apply such a procedure for damage identification, especially for automatic damage detection. Jiang et al. [8] applied multiwavelet packet as a prefilter to enhance the weak multifault features in the narrow frequency bands. Then two ranges were suggested for the amplitude: (0–0.2) of the standard deviation of the original signal for high frequency components and (0.2–0.6) of the standard deviation of the original signal for the low frequency components. As some specific amplitudes were selected (0.04, 0.08, and 0.5) without any justification in this study, it seems that no robust guide is yet available to choose the optimum amplitude based on the wide suggested ranges. Tabrizi et al. [9] applied the wavelet packet decomposition with combination of the EEMD to identify very small faults under various operating conditions. It was concluded that more appropriate amplitude was (0.4–0.6) of the standard deviation for noisy signals and 0.5 for denoised signals. The number of trial was set on 100 for all conditions.

As we have shown, there is no reliable guide for amplitude settling; so far, in this study, a new method (PIEEMD) is proposed to calculate an appropriate and effective noise amplitude for real vibration signals. A number of vibration signals is analysed to verify the proposed algorithm in automatic fault diagnosis based on support vector machine (SVM). Furthermore, as mentioned before, there is no suggestion on the specific number for ensemble trials; this is only declared in three aforementioned studies: the modified new number for simulated signals (70 trials) in [5], a very high number (3000 trials) in [6], and the 100 trials in [3, 9]. In view of this, the appropriate ensemble trial number for real data is investigated as well.

2. EMD Algorithm

The EMD method decomposes a complex signal into a number of IMFs. Decomposition consists of following steps [2].(1)To identify all the local extrema and then connect all the local maxima by an interpolation method. To repeat the procedure for the local minima to produce the lower envelope.(2)To determine the difference between the signal and the mean of upper and lower envelope value to obtain the first component. If it is an IMF, then it would be the first component of . Otherwise, it is treated as the original signal and steps - are repeated. The sifting process can be stopped by any of the predetermined criteria which will be discussed in Section 3.(3)To separate IMF from the original signal to obtain the residue and consider it as the new data and repeat the above described process.(4)To stop the decomposition process when the residue becomes a monotonic function from which no more IMF can be extracted.

3. Ensemble Empirical Mode Decomposition (EEMD)

Decomposition using EEMD consists of the following steps.(a)To add a random white noise signal to the acquired original signal:where Amp is the amplitude of added white noise and is the predetermined number of trial.(b)To decompose the obtained signal () into IMFs using EMD:where represents the th IMF of the th trial, represents the residue of th trial, and is the IMFs number of the th trial.(c)To repeat steps (a) and (b) until the predefined ensemble trial number () (add different random noise signal each time).(d)To calculate the ensemble means of the corresponding IMFs of the decompositions as the final result (): where is the minimum number of IMFs among all the trials.

4. Performance Improved Ensemble Empirical Mode Decomposition (PIEEMD)

As mentioned in Section 1, the added noise must affect the extrema of the original signal so that the intermittency of the components will be removed or decreased as much as possible. However, in the predefined constant amplitude value, the extrema are being affected (and as a consequence decreasing the existed mode mixing) by a random noise, which might not effectively change some extrema.

Instead, an adaptive method (PIEEMD) is proposed and its performance and applicability are evaluated utilizing several real vibration signals. After adding a random white noise, by applying the signal-to-noise ratio (SNR) definition (4), the Amplitude value for each data point of a sample is obtained from (5). Considering an appropriate value for SNR, there would be a confidence that the extrema of the original signal are influenced adequately:where .

In Figure 1, a vibration signal of a roller bearing and a created random noise are shown. A suggested fixed value (0.3) multiplied by standard deviation of the original signal creates a predefined constant value along the whole signal (Figure 2). Thus, affecting the extrema depends on value of random noise at the location of the extrema. Using the proposed algorithm (5), an adaptive value (Figure 2) is generated to preserve the SNR ratio. It means that, for any randomly created noise, the amplitude will be high enough to affect the extrema. Investigating the result of adding noise to the vibration signal shows how the proposed amplitude acts more efficiently on the extrema (Figure 3).

Figure 1: A real collected vibration signal from a roller bearing (a) and a created random white noise (b).
Figure 2: The predefined constant amplitude (0.3 std (signal)) (a) and the proposed amplitude algorithm (PIEEMD with SNR = 10) (b).
Figure 3: The influence of using constant amplitude (0.3 std (signal)) and the PIEEMD algorithm (SNR = 10) on extrema.

5. Support Vector Machine (SVM)

Support vector machine (SVM) is a powerful technique for data classification [10]. Based on available acquired data (training data), the SVM attempts to construct a hyperplane that separates two different classes of samples and orients it to maximize the “Margin” which is the distance from the hyperplane to the closest data points in either class. An example of the optimal hyperplane of two datasets is shown in Figure 4 from [11]. Every time a new element appears, it can be classified according to where it places with respect to the separating hyperplane. The SVM could also be applied in a case of nonlinear classification by mapping the data onto a high dimensional feature space (), where the linear classification is then possible. By applying Kernel function as the inner product of mapping functions it is not necessary to explicitly evaluate mapping in the feature space. Various kernel functions could be used, such as linear, polynomial, or Gaussian RBF (radial basis function). In real world problem, it might be impossible to get a separate line dividing the data and a curved decision boundary might be chosen inevitably. However, ignoring few outlier data points will create smooth boundary. This is handled here by using the slack variable and the error penalty ; it is called soft margin-SVM.

Figure 4: Classification of two classes of data using SVM.

The Margin is defined as [12] And the optimization problem will be as follows [12]: where and are the vector and scalar that are used to define the position of the hyperplane and is measuring the distance between the hyperplane and the samples laying in the wrong side of the hyperplane. Introducing Lagrange multipliers and solving the dual optimization problem, nonlinear decision function will be [11]

6. Methodology

The goal of this study is to evaluate the performance of the proposed amplitude calculation algorithm (PIEEMD) for the various operating conditions of a roller bearing.

The fault diagnosis method is given as follows.(1)To collect vibration signals both for healthy and defective bearings at three different external loads and four shaft speeds.(2)To apply the EEMD with different amplitude of added white noise to decompose the vibration signals into some IMFs. The first IMFs including the most dominant fault information are chosen to extract the feature.(3)To calculate the total energy of the first IMFs:(4)To create a feature vector with the energies of the selected IMFs:(5)To normalize the feature: where .(6)To carry out the training procedure of SVM by utilizing the normalized feature vectors. The 60% of data are used for training and the rest are taken as the test samples.After training the SVM successfully, it would be ready to test samples to identify the different work conditions and fault patterns.

7. Experiment

The bearing dataset (acceleration signals) was collected from an aeronautical specific bearing (SKF NN 76141) under various operating conditions using the test rig (Figure 5) developed and assembled by the Dynamics and Identification Research Group (DIRG) at the Department of Mechanical and Aerospace Engineering of Politecnico di Torino. The signals were acquired by Oros OR38 data acquisition system at 102.4 kHz sampling frequency with a passing band of 40 kHz, taking into account the antialiasing filters. The Kistler triaxial accelerometers (model 8763A500) were used, in which the axes orientation is shown in Figure 5 so that x-, y- and z-axis correspond to the axial, radial, and tangential direction, respectively. In addition to the healthy bearing, two defective bearings were utilized during the test, one with the very small artificial defect severity over one roller (150 microns in diameter) and another with the same fault level on the inner ring. Four different shaft speeds (100, 200, 300, and 400 Hz) and three different external radial loads (1.0, 1.4, and 1.8 kN) were considered to acquire the signals in different operating and controlled conditions, allowing speed, load, and oil temperature control.

Figure 5: DIRG test rig, the axes orientation of the triaxial accelerometers (x, y, and z) and the damaged roller used in the tests.

The original acquired signals were divided into 20 segments including 10000 data points each, to extract required informative feature vectors. Thus, each signal includes 20 segments which create 20 feature vectors as inputs for the SVM. Selecting samples as the training ones includes all the possible random selections to obtain the maximum classification accuracy rate for training.

8. Analysis

Implementing the methodology described in Section 5, the feature vectors for each algorithm, damage location, and signal direction are obtained. The normalized energy of IMFs introduced as an efficient feature vector in fault diagnosis of roller bearing has been adopted just using only first three elements of the feature vectors [13].

In Table 1, it is shown that the smallest defect size (150 microns on a roller) is not recognized using EMD for signals collected through the accelerometer in Y direction (for some working conditions). Although applying the EEMD (with the 100 trails) improves the success rates, there is no correct classification and fault diagnosis for some operating conditions. It seems that amplitudes with 0.3, 0.5, and 0.6 lead to less misclassification (only three operating conditions).

Table 1: The classification rates for different operating conditions (signals collected in direction).

Now, by applying the PIEEMD algorithm, the accuracy of damage detection is investigated. The success rate of defect detection is shown in Table 2 for various preselected SNR values. Obviously, a considerable improved success rates are achieved for some SNR values, especially for SNR = 10 so that there exists only one working condition (speed = 200 Hz and load = 1.8 kN) that the state of bearing is not perfectly identified. However, none of the SNR values leads to a perfect labelling for all conditions. Increasing the SNR (to SNR = 20) increases the success rate for the mentioned operating condition which means that such a signal needs weaker noise to affect the extrema and decrease the mode mixing. Whereas for some signals (such as speed = 100 Hz/load = 1.8 kN and speed = 300 Hz/load = 1.4 kN) the smaller value (SNR = 10) seems to be more appropriate. It means that those signals require some stronger noises.

Table 2: The success rates of damage detection using the proposed algorithm (PIEEMD) in direction.

Exploring reliability of the obtained success rate, the Margin ((6) and Figure 4) of each SVM classification is calculated and presented in Table 3. It is clear from the definition of Margin that higher Margin means more reliable hyperplane and classification. As it can be seen, some Margins are much smaller than others such as 0.6806 (SNR = 20, speed = 300 Hz, and load = 1.4 kN). It means that the possibility of sorting out new misclassified samples (like the result shown in Table 2) is increased. On the other hand, it is reasonable to expect that higher Margins have more reliable results (a correct classification and defect detection for any new investigated sample). The most important conditions are those achieving a perfect classification rate (100%), whereas the calculated Margin is not high enough such as 0.7749 (SNR = 20, speed = 200 Hz, and load = 1.8 kN). There exists the possibility of misclassification for the new samples based on the constructed SVM. It is worth mentioning that there is no determined reliable value for Margin. The Margins calculated are shown in Table 4 for the amplitudes who had better results.

Table 3: The Margin calculated using the proposed algorithm (PIEEMD) in direction.
Table 4: The Margin calculated using EEMD with different amplitudes ( direction).

To test the reliability of the constructed SVM, 20 new samples (10 healthy and 10 damaged samples) for each operating condition were classified with previously constructed SVM (Table 2). The results are proposed in Table 5. Obviously, as it is expected, the new samples are not classified perfectly for previously mentioned low Margin (SNR = 20, speed = 200 Hz, and load = 1.8 kN). Although the previous success rate was 100% and it seemed to be a reliable constructed SVM, its low Margin (in comparison with those showing a perfect damage detection) declares that it may not be a confident SVM (Figure 6). However, for all other conditions, which had higher Margin, the states of the new samples are identified correctly. As it can be seen in Figure 7, all faulty and healthy samples are completely separable.

Table 5: The reliability test of the constructed SVM with 20 new samples ( direction).
Figure 6: The normalized energy of the three first IMFs () with new samples, SNR = 20 (speed = 200 Hz and load = 1.8 kN).
Figure 7: The normalized energy of the three first IMFs () with new samples, SNR = 10 (speed = 200 Hz and load = 1.0 kN).

For the EMD and EEMD (with the predefined constant amplitudes or the presented method (SNR = 10)), the previously constructed SVMs can successfully detect the fault, as it was expected because of their high Margins. The results shown in Table 5 confirm such an expectation.

To investigate in more detail, the collected signals of another defective bearing (small defect on the inner ring) in two directions (Y and Z) are analyzed. The results of classification are shown in Table 6. As all the constructed SVMs have high Margins, they are reliable and lead to perfect success rates for both Y and Z directions, except for one condition (speed = 200 Hz and load = 1.8 kN) in Z-axis which achieved 81.3% success rate.

Table 6: The Margin calculated (using PIEEMD with SNR = 10), for the defective inner ring.

Finally, the appropriate ensemble trial number for damage detection is investigated. The results of fault diagnosis for different trials are presented in Table 7. In the most of the operating conditions, the 30 ensemble trial numbers are enough to achieve a perfect labelling rate, which is far lower than the number used in previous studies (as mentioned in Section 1). It leads substantially to less time consuming and decreasing the computational cost. The higher Margin is not obtained by increasing the number; even in some cases the 30 trials achieve a little higher Margin. For two conditions (speed = 200 Hz, load = 1.4 kN and speed = 300 Hz, load = 1.4 kN), the 100% success in detection is achieved by applying the 100 trails. There exists only one operating condition in which the 100 trials do not achieve perfect success rate (speed = 200 Hz and load = 1.8 kN). As it can be seen in Table 8, implementing the proposed method with 50 trials improves the success rate and the Margin. Increasing the trials to 1000, higher classification accuracy is achieved (from 87.5% to 93.8%), whereas the Margin does not change considerably after the 50 trial numbers. Although with the 70 trials the classification result improves from 81.3% to 87.5%; however, the original EEMD does not obtain higher rates by applying more trial numbers (even up to 1000).

Table 7: The results of using different trial numbers (NR) using PIEEMD (SNR = 10) for direction.
Table 8: The results of using different trial numbers (NR) for the condition: speed = 200 Hz and load = 1.8 kN.

On the contrary, denoising the signals can improve the results to obtain perfect classification results [9]. Using only 30 trials, the same success rate (93.8%), achieved with 1000 trials for the noisy signals, can be obtained (Table 9). To achieve 100% labelling success, the limit of 50 trials is sufficient; however, there is no considerable improving of the Margins after the 50 trial numbers.

Table 9: The results of using different trial numbers (NR) after denoising, for the condition: speed = 200 Hz and load = 1.8 kN.

9. Conclusion

Obviously, in the EEMD technique, there is no robust guide to select an optimal amplitude for the added noise. As the amount of the amplitude is very important, especially for early damage detection (very small defects), in this study a reliable method (PIEEMD) is proposed to calculate the amplitude, instead of using previously suggested amplitudes, a priori predefined and constant. The vibration signals of the various operating conditions are analyzed for three bearings with different states: healthy, faulty (small fault on a roller), and faulty (small defect on inner ring). It is shown that, by applying the proposed amplitude calculating algorithm (especially with SNR = 10), considerable improvement in accuracy of damage detection is obtained, in comparison with predefined constant amplitudes. Exploring reliability of the obtained success rates, the Margin of each SVM classification is calculated and it is confirmed that, for those conditions whose Margin is relatively high, the results are more reliable. For the defective inner ring, the acceleration signals of two radial directions are investigated to achieve more confident results and it is validated that the PIEEMD algorithm looks reliable and can be favorably applied instead of the previous predetermined approach.

It is shown that substantially low trial numbers are required to achieve perfect labelling of samples, which leads to considerably less computational cost. In the most operating conditions only the 30 trial numbers are sufficient to achieve successful damage detection. Using the 100 trial numbers lead to 100% success rate for all working conditions, except for one condition (speed = 200 Hz, load = 1.4 kN) that although increasing the number improves its result, the perfect classification can be achieved only after denoising.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This work has been partially carried out in the framework of the GREAT2020-phase II, project.

References

  1. R. B. Randall and J. Antoni, “Rolling element bearing diagnostics—a tutorial,” Mechanical Systems and Signal Processing, vol. 25, no. 2, pp. 485–520, 2011. View at Publisher · View at Google Scholar · View at Scopus
  2. N. E. Huang, Z. Shen, S. R. Long et al., “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society of London Series A, vol. 454, no. 1971, pp. 903–995, 1998. View at Publisher · View at Google Scholar · View at MathSciNet
  3. Z. Wu and N. E. Huang, “Ensemble empirical mode decomposition: a noise-assisted data analysis method,” Advances in Adaptive Data Analysis, vol. 1, no. 1, pp. 1–41, 2009. View at Publisher · View at Google Scholar · View at Scopus
  4. Y. Lei, J. Lin, Z. He, and M. J. Zuo, “A review on empirical mode decomposition in fault diagnosis of rotating machinery,” Mechanical Systems and Signal Processing, vol. 35, no. 1-2, pp. 108–126, 2013. View at Publisher · View at Google Scholar · View at Scopus
  5. J. Zhang, R. Yan, R. X. Gao, and Z. Feng, “Performance enhancement of ensemble empirical mode decomposition,” Mechanical Systems and Signal Processing, vol. 24, no. 7, pp. 2104–2123, 2010. View at Publisher · View at Google Scholar · View at Scopus
  6. W. Guo and P. W. Tse, “Enhancing the ability of ensemble empirical mode decomposition in machine fault diagnosis,” in Proceedings of the Prognostics and System Health Management Conference, January 2010. View at Publisher · View at Google Scholar · View at Scopus
  7. J. Lin, “Improved ensemble empirical mode decomposition and its applications to gearbox fault signal processing,” International Journal of Computer Science Issues, vol. 9, no. 6, pp. 194–199, 2012. View at Google Scholar
  8. H. Jiang, C. Li, and H. Li, “An improved EEMD with multiwavelet packet for rotating machinery multi-fault diagnosis,” Mechanical Systems and Signal Processing, vol. 36, no. 2, pp. 225–239, 2013. View at Publisher · View at Google Scholar · View at Scopus
  9. A. Tabrizi, L. Garibaldi, A. Fasana, and S. Marchesiello, “Early damage detection of roller bearings using wavelet packet decomposition, ensemble empirical mode decomposition and support vector machine,” Meccanica, vol. 50, no. 3, pp. 865–874, 2015. View at Publisher · View at Google Scholar · View at Scopus
  10. V. N. Vapnik, The Nature of Statistical Learning Theory, Springer, Berlin, Germany, 1995. View at Publisher · View at Google Scholar · View at MathSciNet
  11. B.-S. Yang, T. Han, and W.-W. Hwang, “Fault diagnosis of rotating machinery based on multi-class support vector machines,” Journal of Mechanical Science and Technology, vol. 19, no. 3, pp. 846–859, 2005. View at Publisher · View at Google Scholar · View at Scopus
  12. C. J. C. Burges, “A tutorial on support vector machines for pattern recognition,” Data Mining and Knowledge Discovery, vol. 2, no. 2, pp. 121–167, 1998. View at Publisher · View at Google Scholar · View at Scopus
  13. A. Tabrizi, L. Garibaldi, A. Fasana, and S. Marchesiello, “Influence of stopping criterion for sifting process of empirical mode decomposition (EMD) on roller bearing fault diagnosis,” in Advances in Condition Monitoring of Machinery in Non-Stationary Operations, G. Dalpiaz, R. Rubini, G. D'Elia et al., Eds., Lecture Notes in Mechanical Engineering, pp. 389–398, Springer, Berlin, Germany, 2014. View at Publisher · View at Google Scholar