Emmanuel Ramasso
email: emmanuel.ramasso@femto-st.fr.

Vincent Placet
email: vincent.placet@univ-fcomte.fr

come

Unsupervised consensus clustering of acoustic emission time-series for robust damage sequence estimation in composites

Unsupervised consensus clustering of acoustic emission time-series for robust damage sequence estimation in composites

Emmanuel Ramasso, Vincent Placet, M. Lamine Boubakar

Abstract-This paper suggests a new approach for unsupervised pattern recognition in acoustic emission (AE) time-series issued from composite materials. The originality holds in the development of a clustering ensemble method able to emphasize sudden growths of damages in composites under sollicitations. The method combines multiple partitions issued from different parametrizations, initial conditions and algorithms. A first stage automatically select multifarious subsets of features based on the entropy of sequences of damages detected by clustering. A polygonal representation of the sequences is suggested in order to emphasize the kinetics of fracture events. The second stage allows to estimating the optimal number of clusters necessary to represent the structure of the AE data stream. The data structure is estimated by consensus clustering with boostrap ensembles, which allows to estimating the uncertainty envelopes of each cluster and giving access to an interval of cumulated loading thresholds necessary to activate a particular damage. A qualitative evaluation phase is proposed on simulated datasets to statistically assess and underline both the robustness and accuracy of the proposed clustering fusion method, comparing Kmeans, Gustafson-Kessel algorithm and Hidden Markov Models. An application is then presented for the detection of early signs of failure in high performance carbon fibre-reinforced thermoset matrix composites dedicated to severe operating conditions. Despite the complexity of the configuration (ring-shaped specimens, high emissivity), it is demonstrated that the method emphasizes damage onsets and kinetics (fibre tow breakage, hoop splitting and delamination) within the unevenly-spaced AE timeseries recorded during loading.

Index Terms-acoustic emission, composite materials, clustering fusion, number of cluster estimation, feature selection, uncertainty quantification, structural health monitoring, unevenlyspaced time-series.

I. INTRODUCTION

The health monitoring of critical components in engineering systems is of paramount importance to ensure equipment safety and availability. Those components are subjected to the action of external stimuli which may be recorded by different non-destructive techniques (NDT). This paper is focused on health monitoring based on acoustic emission (AE), a technique particularly well suited for real-time monitoring of the transient elastic waves produced by the sudden redistribution of stress in materials due to local fractures propagating to the surface [START_REF]ASTM Standard Definitions of Terms Relating to Acoustic Emission[END_REF]. The AE technique was used in many applications such as civil structures [START_REF] Behnia | Advanced structural health monitoring of concrete structures with the aid of acoustic emission[END_REF], storage in tanks [START_REF] Kwon | Acoustic emission testing of repaired storage tank[END_REF], pipes corrosion [START_REF] De La Rosa | Higher order statistics and independent component analysis for spectral characterization of acoustic emission signals in steel pipes[END_REF], gears in helicopters, wind turbines and other rotational machines [START_REF] Ruoyu | Rotational machine health monitoring and fault detection using EMD-based acoustic emission feature quantification[END_REF], [START_REF] Qu | Gearbox tooth cut fault diagnostics using acoustic emission and vibration sensors -A comparative study[END_REF], bearings [START_REF] Kilundu | Cyclostationarity of acoustic emissions (AE) for monitoring bearing defects[END_REF], cutting tools [START_REF] Martins | Tool condition monitoring of single-point dresser using acoustic emission and neural networks models[END_REF], electronics [START_REF] Ukil | Monitoring of HV generator circuit breaker contact ablation based on acoustic emission[END_REF] and medical diagnosis [START_REF] Mascaro | Exploratory study of a non-invasive method based on acoustic emission for assessing the dynamic integrity of knee joints[END_REF]. This paper is particularly focused on the detection of damages in composite structures. Structural health monitoring (SHM) based on AE is of paramount importance for the transportation and aerospace industries with the aim to reduce energy comsumption while improving equipment performance.

AE-based SHM involves to face one important difficulty that is to discriminate, precisely and with robustness, AE signals generated by a specific damage from other AE sources. In composite materials, three main families are generally considered, namely matrix cracking, fibre breakages and interface failures [START_REF] De Groot | Real-time frequency determination of acoustic emission for different fracture mechanisms in carbon/epoxy composites[END_REF], [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]. Among these families, sub-families can also be distinguished. External sources are also generally present generating mixed transients due to actuating systems and to the testing environment that can generate electromagnetic and mechanical noise, rubbing and friction [START_REF] Harris | Fatigue and accumulation of damage in reinforced plastics[END_REF], [START_REF] Barat | Detection of AE signals against background friction[END_REF], [START_REF] Cheng | Internal friction behavior of unidirectional carbon/carbon composites after different fatigue cycles[END_REF]. The discrimination of AE sources based on AE signals is generally performed by data mining techniques. Usual approaches for mining AE data are either based on supervised [START_REF] Otsu | Pattern recognition analysis of acoustic emission from unidirectional carbon fiber-epoxy composites by using autoregressive modeling[END_REF], [START_REF] Godin | Clustering of acoustic emission signals collected during tensile tests on unidirectional glass/polyester composite using supervised and unsupervised classifiers[END_REF], unsupervised [START_REF] Murthy | Application of pattern recognition concepts to acoustic emission signals analysis[END_REF], [START_REF] Gutkin | On acoustic emission for failure investigation in CFRP: Pattern recognition and peak frequency analyses[END_REF], [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF], [START_REF] Doan | An unsupervised pattern recognition approach for ae data originating from fatigue tests on polymer-composite materials[END_REF] or partially-supervised [START_REF] Ramasso | Health assessment of composite structures in unconstrained environments using partially supervised pattern recognition tools[END_REF] classification tools according to the amount of prior information available.

In non-trivial experimental configurations similar to industrial in-service conditions, the monitoring of a particular damage mode using AE is particularly challenging. This is the case with highly emissive composite materials that produce a massive amount of AE signals under sollicitations [START_REF] Doan | An unsupervised pattern recognition approach for ae data originating from fatigue tests on polymer-composite materials[END_REF]. Non-planar geometry such as tubular structures [START_REF] Capela | Fatigue strength of tubular carbon fibre composites under bending/torsion loading[END_REF] used in energy storage and underwater applications have both complex stress state and damage process involving complicated AE waveforms. In-service conditions also involve non-stationary noise which is superimposed on AE signals leading complex AE signals [START_REF] Kharrat | A signal processing method for hits detection and separation in high AE activity systems: Application to composite materials under fatigue tests[END_REF]. In addition, cumulated damage in composites and the associated discontinuities, also induces a distortion of transients propagating in the material [START_REF] Sause | Acoustic emission signal propagation in damaged composite structures[END_REF], [START_REF] Kharrat | Influence de l'endommagement accumulé sur la signature acoustique au cours d'un essai de traction sur un composite[END_REF]. Therefore, even if the prototypes of damages are initially compact, their characteristics can evolve in the feature space and their shape may change as a function of damage level. For such non-trivial configurations, unsupervised data mining techniques are more suitable.

From a clustering point of view, the first difficulty concerns the scattering in clustering results that can be induced by algorithm parametrization. Indeed, only a single set of parameters are generally used in existing algorithms to estimate the clusters from complex AE datasets. The second difficulty concerns the scattering of AE features within clusters due to the aforementioned causes.

The algorithms developed in the literature for unsupervised AE-based SHM generally do not take into account the kinetics of damages. It relies on the application of standard algorithms (Kmeans, Fuzzy-C-means, Gaussian Mixture Models, Self-Organising Maps, Gustafson-Kessel algorithm) [START_REF] Sawan | Unsupervised learning for classification of acoustic emission events from tensile and bending experiments with open-hole carbon fiber composite samples[END_REF] that only allow point-wise evaluation of clusters without explicit correlation between points in space or in time. The aim of those algorithms is generally stated as finding the most compact and separated clusters [START_REF] Kostopoulos | Fracture behavior and damage mechanisms identification of SiC/glass ceramic composites using AE monitoring[END_REF], [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF], [START_REF] Pomponi | A real-time approach to acoustic emission clustering[END_REF].

Although most of damages in composite materials occur sequentially or in cascade [START_REF] Pollock | Acoustic Emission and Critical Phenomena: From Structural Mechanics to Geophysics[END_REF], only a few clustering methods proposed in the literature were focused on the temporal evolution of clusters estimated from AE data. In [START_REF] Kostopoulos | Fracture behavior and damage mechanisms identification of SiC/glass ceramic composites using AE monitoring[END_REF], the authors were specifically interested in the triggering of specific damages and the evolution in time of the different damage mechanisms in composites. For that, they studied distance measures to discover compact clusters. In [START_REF] Dzenis | Analysis of microdamage evolution histories in composites[END_REF], the shape of the frequency spectrum was used to create libraries of damages. Those libraries allowed to classifying AE data, and the evolution in time of clusters was correlated to the cumulated number of AE counts and to the distribution of amplitudes. The Kmeans algorithm was used in [START_REF] Silversides | Acoustic emission monitoring of interlaminar delamination onset in carbon fibre composites[END_REF] to build a noise model in order to discriminate between both fatigue signal noise and acoustic emission signals emitted from crack initiation and crack growth during delamination. Evolving clustering methods were proposed in [START_REF] Pomponi | A real-time approach to acoustic emission clustering[END_REF], [START_REF] Placet | Online segmentation of acoustic emission data streams for detection of damages in composites structures in unconstrained environments[END_REF], [START_REF] Doan | An unsupervised pattern recognition approach for ae data originating from fatigue tests on polymer-composite materials[END_REF] where new clusters can be added dynamically according to the characteristics of AE hits observed during quasi-static and fatigue loading. Statistical modelling of AE time-series based on Hidden Markov Models (HMM) was suggested to estimate the temporal evolution of clusters [START_REF] Zhou | Damage classification for structural health monitoring using time-frequency feature extraction and continuous hidden Markov models[END_REF]. The problem is that AE time-series are unevenly-spaced, as encountered in geoscientific measurements [START_REF] Rehfeld | Comparison of correlation analysis techniques for irregularly sampled time series[END_REF], meaning that the spacing between AE hits is not constant. Therefore specific algorithms may be required for the post-treatment of the partitions into relevant clusters [START_REF] Ramasso | Health assessment of composite structures in unconstrained environments using partially supervised pattern recognition tools[END_REF].

This paper suggests a new approach for the clustering of AE time-series with an attempt to estimate the sequence of damages without assumption about the compacity of clusters. For that, a method is proposed which has the following characteristics: 1) Optimally select the number of clusters without using additional distance measures for partition evaluation but using a probabilistic formulation; 2) Combine multiple clusterings obtained from multifarious subsets instead of a single subset of features; 3) Represent the uncertainty around the estimation of both the damage onsets and the loading thresholds necessary to activate those damages. Robustnessto-parameterization of the proposed method is attained by clustering fusion, also called clustering ensemble or consensus clustering in the literature [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF]. A particular attention is paid to the Gustafson-Kessel clustering algorithm [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] to cope with data scattering which is of particular importance for SHM applications [START_REF] Farrar | Structural Health Monitoring: A Machine Learning Perspective[END_REF].

Section II is dedicated to the presentation of the method. Section III is focused on illustrations of each step of the method and to a comparison with other methods. The application on a ring-shaped composite specimen made of carbon fibres and thermoset resin is then described in Section IV.

II. ENSEMBLE METHOD FOR ACOUSTIC EMISSION

CLUSTERING

The method is decomposed in two phases (Figure 1), each illustrated in the section dedicated to results. The first phase (Section II-A) is dedicated to the selection of multifarious subsets; The second phase (Section II-B) to the optimal selection of the number of clusters and to uncertainty quantification about the sequence of fracture events. The method considers a set of L parameterizations (such as subsets of features or algorithms) S 1 , S 2 . . . S L generating hard partitions with K clusters.

1) Cumulated activity of AE sources (Step 1A): Given a hard partition P t ∈ {1 . . . K}, t = 1 . . . T , the logarithm of the Cumulated Sum of a Cluster Appearance of a cluster j over time (or loading), CSCA jt , is computed to emphasize the dynamics of the activity of an AE source:

L jt = log t t =1 I P t =j , t = 1 . . . T (1)
where the indicator function I Pt=j ≡ 1 if cluster j is activated at time t. The set of K CSCA curves is called damage profile.

For a given K, the clusters are sorted with report to the final value L j,T for each subset. Sorting allows fast co-association without computing a T ×T co-association matrix [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF] that can be huge according to the emissivity of the material considered. By sorting, two clusters C i in a partition P a and C j in a partition P b are co-associated if both clusters have the same rank.

The parameterizations that lead to clusters with different kinetics are prefered. For that, the area between two consecutive curves L j,1:T and L j+1,1:T , j = 1 . . . K is used to discriminate clusters in terms of kinetics. The couples (t 1 , L jt1 ), (t 2 , L jt2 ), . . . (t T , L jT ), (t T , 0), (t T -1 , 0), . . . (t 2 , 0), (t 1 , 0), (t 1 , L jt1 ) are considered as ordered pairs of vertices listed counterclockwise and representing a polygon (the last point allows to closing it) [START_REF] Rosen | Handbook of discrete and computational geometry[END_REF], [START_REF] Ramasso | Investigating computational geometry for failure prognostics[END_REF]. The polygonal representation P(j) of the j-th cluster has a measurable area that can be computed efficiently using the Shoelace formula [START_REF] Braden | The surveyor's area formula[END_REF]. The area A(j, j + 1) between two curves related to two clusters j and j + 1 is then given by the difference between the area of both polygons:

A(j, j + 1) = P(j) -P(j + 1), j = 1 . . . K -1 (2)
The aforementioned procedure is applied for each parameterization S 1 . . . S l , . . . S L .

2) Outlier removal (Step 1B): An outlier removal procedure is developed to decrease the number of potential parameterizations used further in the fusion process. For that, the trend of successive areas between clusters (sorted by descending order according to the number of AE signals) is estimated by a polynomial of degree 2 that best fits the areas

A(1, 2), A(2, 3), . . . A(K -1, K) in a least-squares sense: A(j, j + 1) ≈ p 2 • j 2 + p 1 • j + p 0 , j = 1 . . . K -1 (3)
Coefficient p 2 represents the shape of the parabola and gives the trend of the areas (increasing or decreasing). The outlier removal thus consists in keeping only the parameterizations belonging to the set of statistically relevant ones, S relevant , estimated by the range of percentiles (IPR) around the median:

S relevant = {S l | p l 2 ∈ IP R}, l = 1 . . . L (4)
with

IP R = Q 0.5-β (p 1:L 2 ), Q 0.5+β (p 1:L 2 ) , 0 < β < 0.5 (5) where Q n (p 1:L 1 )
represents the quantile for probability n of the values of the parabolas' shape (p 2 ) for all parameterizations (n = 0.5 corresponds to a median). In experiments, β = 0.25 leading to the interquartile range.

3) Entropy-based parameterization ordering (Step 1C): The redundancy of clusters' kinetics obtained with the l-th parameterization (l = 1 . . . L, L = |S relevant |) is estimated by the entropy of the areas:

E l = - K-1 k=1 A l * (k, k + 1) log 2 A l * (k, k + 1), l = 1 . . . L (6)
where A l * is the normalised area given by:

A l * (k, k + 1) = A l (k, k + 1) K-1 k =1 A l (k , k + 1) , l = 1 . . . L (7)
Following the Principle of Maximum Entropy, the least committed probability distributions over possible clusters are selected by considering those above the 90-th percentile (Q 0.90 (E 1:L )):

S final = {S l ∈ S relevant , l = 1 . . . L | E l ≥ Q 0.90 (E 1:L )} (8)
The threshold (here the 90-th percentile) can be adjusted to obtain about 20 -50 partitions which can be practically sufficient for fusion [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF].

4) Partition fusion (Step 1D):

The clusters have a similar semantic on all partitions due to sorting and co-association (Section II-A1). It follows that the majority voting rule can be applied that assigns to an AE signal the cluster receiving the highest number of votes. Given the partitions P l t ∈ {1, . . . K}, t = 1 . . . T, l = 1 . . . L, L = |S final |, the number of votes for a cluster j is given by:

N j = L l=1 I P l t =j , l = 1 . . . L (9)
The final decision on the cluster membership of a data point x t (feature vector attached to an AE signal) is given by:

P * t = argmax j=1...K N j (10)
B. Optimal choice of K and quantification of uncertainty (Phase 2)

Each of the L parameterizations selected in the previous section (S final ) leads to one partition P l,K , and the fusion to one additional combined partition P * K , each made of K clusters.

1) Optimal number of clusters (Step 2A): It is assumed that the selected subsets of features obtained for K are relevant whatever the number of clusters. Therefore, the "optimal" number of clusters can be selected in three steps:

1) Vary the number of clusters k from k min to k max and compute the partitions P l,k , l = 1 . . . L ; 2) Apply the fusion process and compute P * k ; 3) Select the "best" partition P * k k k * and deduce the value of k * . The optimality criterion proposed in [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF], [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF] is used to select the optimal value of k and to quantify the robustness of the fusion process with report to small changes in the initial set of partitions. The method considers that the clustering ensemble

P = {P l,k , k ∈ [k min , k max ] , l = 1 . . . L} is disturbed by a bootstrap. A set of B bootstrap clustering ensembles P B = {P b1 , P b2 . . . P b B }
is created by sampling with replacement from P. The corresponding set of combined data partitions obtained by the fusion process is denoted as P * B = {P * b1 , P * b2 . . . P * b B }. Then, the consistency of a kcluster combined partition P * k k k B with P B can be estimated by the median value of the Normalized Mutual Information (NMI):

mN M I P * k k k B , P B = B median i=1 N M I P * k k k b i , P bi , (11)
and the NMI between two partitions P a and P b is given by [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF]:

N M I(P a , P b ) = -2 ka i=1 k b j=1 n ab ij log n ab ij •T n a i •n b j ka i=1 n a i log n a i T + k b j=1 n b j log n b j T , (12)
that is different from [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF] in the normalisation. Coefficient n ab ij is the number of shared patterns between clusters C a i ∈ P a and C b j ∈ P b , n a i (resp. n b j ) represents the number of patterns in cluster C a i ∈ P a (resp. C b j ∈ P b ). The NMI is bounded such that 0 ≤ N M I(., .) ≤ 1.

By maximising the median NMI over the number of clusters as follows:

k * = argmax k k k mN M I P * k k k B , P B , (13)
the number of clusters k * can be selected, yielding a "consensus partition" that shares the most information with the clustering ensemble:

P * ≡ P * k k k * B . ( 14
)
The robustness to small variations in the clustering ensemble can be quantified by the interquartile range of NMI values:

IQR k = Q 0.75 N M I P * k k k B , P B -Q 0.25 N M I P * k k k B , P B (15
) where Q 0.25 and Q 0.75 are the first and third quartile respectively of the NMI values.

It is important to notice that the average NMI (in place of the median in Eq. 11) was used in [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF] for the optimal choice of the clusters, whereas the variance (in place of IQR in Eq. 15) was prefered in [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF]. In the latter publication, the authors empirically shown, using some simulated datasets, that the average NMI was biased towards the average number of clusters in the clustering ensemble. This behavior was not observed (in practice) in the present publication, therefore the method proposed in [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF] was chosen.

Compared to usual criterion used in pattern recognition of AE signals such as validity indices, the NMI does not require additional distance measures over features and clusters, enhancing its applicability to various datasets and using different clustering techniques.

2) Uncertainty representation and quantification (Step 2B): The uncertainty of the assignment of clusters applied to pattern discovery in AE signals was investigated recently in [START_REF] Sause | Quantification of the uncertainty of pattern recognition approaches applied to acoustic emission signals[END_REF] using a supervised method based on a ground truth supplied by the end-user. The proposed clustering consensus method makes it possible to represent the uncertainty on clusters in an unsupervised manner. The first way is to estimate the envelope around each cluster represented in the log-CSCA space:

E jt = inf l=1...L L l jt , sup l=1...L L l jt ( 16
)
By superimposing the CSCA of the j-cluster in the final partition P * onto this envelope, it is possible to evaluate the uncertainty of the kinetics. The envelope also gives access to an estimation of an interval ICLT j of the Cumulated Loading Threshold (CLT) necessary to activate cluster j:

ICLT j = CLT j,1:L ± std (CLT j,1:L ) (17)
where CLT j,1:L = {CLT j,l : l = 1 . . . |S final |} is the set of damage onsets of cluster j estimated for each parameterization, CLT j,1:L is the average and std the standard deviation of CLT values given all parameterizations.

The second way to represent uncertainty is to consider accumulated evidence on clusters in a reduced feature space such as the amplitude (in dB) versus duration of AE hits (in µs) [START_REF] Lee | Acoustic emission charaterization of single and dual fiber reinforced metal matrix composites[END_REF], or the second partial power (PP2, in %) versus weighted peak frequency (WPF, in kHz) [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]. For a twodimensional reduced feature space, F 1 (t) versus F 2 (t), t = 1 . . . T , we can proceed as follows: 1) Apply a quantization on both features (for instance by considering fixed intervals with U and V values for each feature) yielding F q 1 (u) and F q 2 (v), u = 1 . . . U, v = 1 . . . V ; 2) For the j-th cluster, count the number of times a region u × v contains a point assigned to cluster j, leading to a map defined by:

M j (u, v) = T t=1 I C1∩C2∩C3 C 1 : x t ∈ F q 1 (u), C 2 : x t ∈ F q 2 (v), C 3 : P * t = j (18)

III. RESULTS ON SIMULATED DATA

The original dataset used subsequently, with the kind permission of M.G.R. Sause, was first published in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF].

A. Damage sequence generation

A data sequence generator is proposed to obtain clusters and AE data constrained by some kinetics. It allows to getting a ground truth used for quantitative evaluation. It is mathematically represented by a three states probabilistic generative model, where each state is supposed to represent one of the three following failure modes: Matrix cracking (s 1 ), interface failure (s 2 ) and fibre breakages (s 3 ). The probability mass assigned to a failure mode at a given time is given up to a constant by:

p(s i ) ≈ 1 1 + exp (-a i • (t -b i )) (19)
where a i and b i determine the shape and the position of the sigmoid functions, normalised to get probability masses. Each sequence is made of 100 data points generated using the following parameters:

a i = 0.4; i = 1, 2, 3; b 1 = 5; b 2 = 25; b 3 = 55; t = 1, 2 . . . 100
For the sake of simplicity, the points are assumed equally spaced (which is not the case for the real data presented in the next Section). The probability mass function is represented in Figure 2. Given each state s i , the data x = [x 1 ; x 2 ; . . . ; x t ; . . . ; x T ], with x t ∈ n a n-dimensional feature vector, are assumed to follow a multivariate Gaussian distribution so that

p(x t |s i ) ∼ N (µ i , Σ i ) (20)
The parameters µ i , Σ i were estimated using the results presented in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF] where the authors studied the degradation of composite plates under four-point bending. The three main failure modes were detected based on AE signals and the Kmeans algorithm. The dataset called B with the 85-th combination (see [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF] for details) was used and the partition was considered as a ground truth so that µ i and Σ i were estimated using the sample mean and covariances in cluster Fig. 2: Probability mass assigned to each failure mode with report to time. The damage level is gradually increasing yielding similar probability of occurrence of all failure modes (for instance matrix cracking, interface failures or fibre breakages).

i. From AE signals obtained in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF], the following n = 14 features were used: the partial power 1 to 6 (PP1 to PP6, in %), representing the proportion of the FFT power spectrum falling in frequency bands [0, 150] kHz, [150, 300] kHz, [300, 450] kHz, [450, 600] kHz, [600, 900] kHz, [900, 1200] kHz, the average frequency, the reverberation frequency, the initiation frequency, the frequency centroid, the peak frequency, the weighted peak frequency, the absolute energy and the amplitude.

A total amount of 100 datasets were generated by using this generative model. For each dataset, Eq. 19 is first applied in order to obtain the state sequence; Then, given the states, the data are generated by sampling using Eq. 20. The generated state sequence is kept in memory to be used as a ground truth (during comparison). Figure 3 represents one dataset in the feature space made of the weighted peak frequency (WPF) and second partial power (PP2). It shows a good agreement between simulated and real data [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]. The true points (a total amount of 1408) are also shown as well as the position of the centers representing the prototypes of failure modes. Weighted Peak Frequency, in kHz; PP2: Second Partial Power, in %, both defined in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]). The real data obtained in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF] are depicted using small gray dots. The figure shows a good agreement between simulated and real data.

B. Illustration of each step of the algorithm

A state sequence (Eq. 19) drawn randomly using the generative model is illustrated in Figure 4(a) using the CSCA representation as described by Eq. 1 (continuous bold lines). One specificity of this representation of clusters concerns the fact that each assignment to a cluster is represented by a step. Multiple consecutive assignments to a cluster may thus represent a sudden growth (avalanche or cascade) of damages which can be observed on composites under sollicitations [START_REF] Pollock | Acoustic Emission and Critical Phenomena: From Structural Mechanics to Geophysics[END_REF]. The aforementioned sequence is accompanied by a 14-dimensional dataset drawn from Eq. 20. All subsets made of 4 and 5 features (3003 subsets representing parameterizations) were considered for clustering using the Gustafson-Kessel algorithm [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF]. The latter was initialised with a Gaussian Mixture Model with 3 components. The fuzziness parameter was set to 1.5 and the j-th cluster volume to ρ j = |Σ j | 1/n as suggested in [START_REF] Kaymak | Fuzzy clustering with volume prototypes and adaptive cluster merging[END_REF] where Σ j is the covariance matrix estimated by the GK algorithm. Each resulting partition was transformed into a CSCA representation (Eq. 1) and the areas were computed (Eq. 2). The polygonal representation and the associated areas of this state sequence is depicted in Figure 4(b). The CSCA obtained for different subsets of features are also represented (grayscale dotted lines) in Figure 4(a). For each subset, the trend of areas between consecutive clusters was estimated by fitting a second degree polynomial (Eq. 3). The 3003 resulting coefficients p 2 are plotted in Figure 5 for all subsets (points) and for the subsets selected by the interquartile range proposed in Eq. 5 with β = 0.25 (circles). Here, only positive values were selected by the algorithm meaning that parabolas mainly opened upward. It means that the areas globally decreased as illustrated in Figure 4(b). Figure 5 pictorially represents the coefficients p 2 selected after sorting and entropy maximization as described in Eq. 6, 7 and 8 (square markers). Grayscale lines in Figure 4(a) represent the sequences obtained by the related subsets. The 90-th percentile (Eq. 8) was here equal to 0.99.

The sequence resulting from the fusion process is illustrated in Figure 4(a). The related partition was compared to the ground truth (continuous bold lines) using the corrected-forchance version of the Rand index called Adjusted Rand Index (ARI) [START_REF] Vinh | Information theoretic measures for clustering comparison: Is a correction for chance necessary?[END_REF] (equals 0 on average for a random partition, and 1 when comparing two identical partitions) yielding 89% for this sequence. The errors (11%) are mainly due to some confusions between state 2 and 3 in [START_REF] Kharrat | Influence de l'endommagement accumulé sur la signature acoustique au cours d'un essai de traction sur un composite[END_REF]55]-th time unit.

Figure 6(a) shows the median of NMI values (Section II-B1) for 100 simulated datasets with report to several number of clusters, using B = 50 boostrap ensembles and 25 subsets of features for each ensemble (drawn randomly from the set S final ). The optimal value (highest median NMI) generally equals 3 which corresponds to the ground truth and the results are robust to variations on features over boostrap ensembles. The performance is evaluated by the number of times the optimal number of clusters is correct, yielding a classification rate of 97% (over 100 datasets). Figure 6(b) depicts the evolution of the IQR of NMI values to estimate the number of clusters yielding 53%. These figures show that the median of NMI is better than IQR to estimate the number of clusters. Moreover, both average and variance proposed in [START_REF] Strehl | Cluster ensembles -a knowledge reuse framework for combining multiple partitions[END_REF] and [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF] led to lower classification rates with 93% and 26% respectively.

The uncertainty is illustrated using the weighted peak frequency feature (WPF, in kHz) and the second partial power (PP2, in %) which were discretized by considering fixed interval of size 10 kHz and 5% respectively. The contours of clusters are represented in Figure 7 by accumulating evidence over a) all subsets of features selected before fusion in phase 1 (left), and b) after fusion with number of clusters optimised (a) Median of NMI values for all datasets (to be maximised). in phase 2 (right). A good agreement can be observed between both the result of the fusion and the model used for generating the data (Fig. 3).

C. Quantitative evaluation

The 100 datasets used in the previous examples are used for comparison purposes. For each dataset, the partitions estimated by three different algorithms were compared to the ground truth using the ARI:

1) The "voting algorithm" proposed in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF] focusing on the selection of a single subset of features on the number of clusters for AE signals ;

2) The proposed method using the Gustafson-Kessel (GK) algorithm for clustering [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF] ; 3) The proposed method using Hidden Markov Models (HMM) for clustering [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]. In the following tests, the main goal was to evaluate the algorithms in subset selection so that K was set to 3 clusters. All subsets made of 4 and 5 features (3003 subsets) were considered. The GK was parameterized as explained previously and the HMM was tuned with 3 states, 1 Gaussian per state and 10 Kmeans-based initialisations [START_REF] Murphy | The Bayes net toolbox for Matlab[END_REF]. Several initialisations are used to quantify the amount of variance in the performance.

Concerning the method proposed in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF], all combinations of 5 features were considered and the Kmeans was applied (using 10 different initialisations and keeping the one minimising the sum of the squared-distances over data points). Seven Fig. 7: Uncertainty on clusters in WPF/PP2 feature space before and after fusion. The fusion is in agreement with the model (Fig. 3).

clustering performance measures were then calculated [START_REF] Nanda | Clustering indian stock market data for portfolio management[END_REF]: the Partition Index, the Separation index, the Xie and Beni's index, the Davies and Bouldin index, the Calinski-Harabasz index, the Krzanowski-Lai index and the Hartigan index. A score was then assigned to each partition for each index as proposed in [START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]: The best partition received 25 points, the second one 24 points and the third one 23 points. The final score was obtained by summing the scores over all indices, the ideal partition was thus expected to receive 7 × 25 points.

Figure 8 illustrates the performances over all subsets using a box plot representation for the three algorithms. On average, the best performance (96%) was obtained for the proposed method using the GK algorithm. When using the HMM, the performance in terms of precision was lower (90%) as well as for the robustness (interquartile range equals 22% for HMM against 9% for GK). The Kmeans-based method led to 12%. The important difference between the two methods (voting based on Kmeans with single feature against multiple clusterings) in terms of performance holds in the fact that the Kmeans optimises the compacity and the circularity of clusters [START_REF] Wu | Top 10 algorithms in data mining[END_REF] which is not relevant for this dataset. On the opposite, the proposed method looks for sequence of clusters. The use of the GK algorithm allows to coping with AE data scattering, where clusters can take the form of ellipsoids, through the use of a Mahalanobis-like distance [START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF]. In addition, the Mahalanobis distance is known to be robust to outliers [START_REF] Peel | Robust mixture modelling using the tdistribution[END_REF]. The proposed method is thus more appropriate when the dynamics follow a sequence of states, when the distribution of features in each state is not globular-shaped and when the number of elements in each cluster is different. It is important to point out that some individual sequences obtained by some subsets of features led to erroneous sequence using the GK clustering algorithm as shown in Figure 4 and in Figures 7(a), 7(c) and 7(e). However, those errors were compensated by the fusion process that draws benefits from multiple clusterings (Figure 7 Fig. 8: Performance assessement: Box plot of the adjusted rand index (ARI) that has to be maximised. The GK coupled with the proposed method for feature selection and estimation of the number of clusters yields the best results in terms of precision (96%) and variation (9%).

IV. EVALUATING THE KINETICS OF DAMAGES ON A RING-SHAPED CFRP COMPOSITE DURING TENSILE TEST Data and results collected on ring-shaped CFRP composite structures tested under quasi-static loading up to failure (total ruine of the material) were investigated and analyzed in this section. Such structures are used in applications related to transportation systems (high-speed rotors, flywheels) under severe operating conditions (involving high and multiaxial mechanical stresses combined with relatively high temperature). The rings were obtained by cutting 6-layers filament wound carbon tubes with lay-up configuration [(90 o ) 2 /±45 o /(90 o ) 2 ]. The split-disk test specimens were machined with respect to the ASTM D2290 standard geometry and dimensions. Each specimen had a width of 23 mm and involved two sections of reduced area with a width of 14 mm. These notches were located 180 o apart from each other with internal diameter equal to 60 mm. Tensile tests were performed using an Instron 8501 testing machine (100 kN) where the specimens were subjected to a displacement control with a 5 mm.min-1 moving crosshead speed.

The transient elastic waves were recorded during test at the material surface using a multi-channels data acquisition system from EPA (Euro Physical Acoustics) corporation (MISTRAS Group). Two miniature piezoelectric sensors (µ-80) were used, characterized by a range of resonance of 250 -325kHz, preamplifiers with a gain of 40 dB, a 20-1000kHz filter, a PCI card with a sampling rate of 2 MHz and the AEWin software for real-time feature extraction with an amplitude threshold at 45 dB. The calibration of the system was performed after installation of the transducers on the specimen and before each test using a pencil-lead break (PLB) procedure. The features defined in [START_REF]PCI-2 based AE system user's manual[END_REF] were used: Average frequency (A freq , in kHz), reverberation frequency (R freq , in kHz), initiation frequency (I freq , in kHz), absolute energy (in aJ), PAC energy (in µV.s), amplitude (A dB in dB), root mean square (mV), Average Signal Level (ASL, in dB), signal strength (in pV.s) and counts-topeak.

For the sake of simplicity, the results on only one specimen are presented in details below. Similar results were observed on multiple composite materials with different shapes (plate, rings) and configuration of layers, and with different matrices (thermoset, thermoplastic) and fibres (carbon and natural).

The number of AE data recorded until failure was 41606. The proposed algorithm was run on all subsets of 4 features (for a total of around 210 subsets) considering 2, 3 . . . 7 clusters. The optimal number of clusters was found to be equal to 4 using the NMI-based criterion.

Figure 9(a) represents the sequences of those four clusters in the CSCA space. Both the individual partitions selected by the algorithm (in gray level) and the result of the fusion (bold lines) are depicted with report to loading. The cumulated energy released during loading is superimposed. From the sequences, it can be observed that:

• The two first clusters started at the beginning of loading. They were the most frequent and characterised by amplitudes around 50 dB and 70 dB (Figure 10(a)). • The third cluster occurred after some occurrences of the two first clusters around 1 kN. This cluster was characterised by relatively high amplitudes (80 -95 dB, Fig. 10(a)). Until 6 kN, it was made of a few appearances represented by a few steps in the CSCA curve. Around 6 kN, it presented a slope change corresponding to a sudden multiple assignments of AE hits to this cluster. This change was synchronised with a slope change observed on the fourth cluster. • The fourth cluster concerned AE hits with the highest amplitudes (above 90 dB, Fig. 10(a)). From 5.5 kN and until failure, theses steps were regularly spaced (with report to loading) meaning that a few AE hits were assigned to this cluster. Some cascades (multiple assignments) occurred around 8.3 kN, 9.7 kN and 13.6 kN. This behavior (order of occurrence and cascades) was observed on different types of composites, with particular kinetics according to the shape, the lay-up configuration and the type of matrix and fibres.

Figure 10(b) shows the number of times each feature was considered in the fusion process. The ASL and the RMS are the most frequent closely followed by the PAC energy, the reverberation frequency, the signal strength and the absolute energy. It can be seen that all features were selected at least once by the algorithm which seems to indicate that all features, represented by multiple subsets, may be necessary to represent the sequence of damages. During loading, the snapshots obtained from the infrared camera (Figures 9(b) and 11) allowed to showing a link between the evolution of the fourth cluster and high releases of heat which are typical of fibre tow breakages. This damage family was observed both during loading by optical observations (Fig. 12) and after the ruine by fractography (Fig. 13). The lay-up configuration implies ply delamination and hoop splitting in 90 o layers that were observed by the optical camera (Fig. 12). The third cluster is likely to be associated to these damages related to interface failures. The second cluster may be related to matrix cracking and also to fretting between the specimen and the clamps of the testing machine (observed post-mortem). The first cluster is attributed to electromechanical noise characterised by AE hits with short duration and low amplitude. (b) Infrared images (front view) recorded at different loading levels, manually annotated using the fuzzy-select GIMP tool in order to enhance the position of highly energetic events (appearing in white) and corresponding to particular activation of some clusters: i) 1.9 kN, ii) 3.7 kN, iii) 8.8 kN, vi) 11.4 kN, v) 12.4 kN, vi) 13.6 kN.

Fig. 9: Sequences of clusters and infrared images obtained for the real composite.

(a) Amplitude range in each cluster. It is important to point out that the proposed method uses multifarious subsets of features with the aim to distinguish between families of damages (clusters) characterised by different kinetics. The chronology is obtained by assigning AE hits into clusters characterised by an order of occurrence in time. Therefore, those families of damages may be composed of multiple elementary damages such as matrix cracking, delamination and fibre breakages. It is in accordance with [START_REF] Ono | Acoustic emission behavior of flawed unidirectional carbon fiber-epoxy composites[END_REF] who showed that acoustic emission signals due to those elementary damages may have overlapping characteristics. Those overlaps may be accentuated by the degradation of the composites during the sollicitations since the Lamb-wave modes can interact with the discontinuities created during the damaging process or due to manufacturing [START_REF] Sause | Acoustic emission signal propagation in damaged composite structures[END_REF].

V. CONCLUSION Unsupervised pattern recognition in AE time-series issued from composite materials was tackled by the use of multiple clusterings. An automatic feature selection was proposed coupled with an optimisation of the number of clusters. Accuracy and robustness were quantified on simulated datasets by using the Gustafson-Kessel algorithm in comparison to two other techniques. The tests demonstrated the importance of using different kinds of features which depict complementarity. Features in both the frequency and the time domain were selected illustrating that the "diversity" brought by those subsets can be beneficial for clustering acoustic emission time series. It is also shown that the methodology is general enough to be used with different clustering methods (Kmeans, Gustafson-Kessel, Hidden Markov Models).

The damage profile estimated by the proposed method also gives access to the main thresholds of damage onsets and kinetics with a quantification of the uncertainties. As illustrated for a real composite, when complex damages occur involving different sources (matrix cracking, fibre-matrix debonding and fibre breakage), the sensitivity of the detection with report to the features can be expectedly high. The proposed method based on clustering fusion allowed to coping with this sensitivity while capturing damage kinetics and onsets.

The coupling of the proposed method with AE source numerical models is under study.

Fig. 1 :

 1 Fig. 1: Plot chart.

Fig. 3 :

 3 Fig.3: Data simulated using the proposed model (WPF: Weighted Peak Frequency, in kHz; PP2: Second Partial Power, in %, both defined in[START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF]). The real data obtained in[START_REF] Sause | Pattern recognition approach to identify natural clusters of acoustic emission signals[END_REF] are depicted using small gray dots. The figure shows a good agreement between simulated and real data.

Fig. 4 :

 4 Fig. 4: State sequence drawn from the model and the estimated sequences represented in the CSCA space. The continuous lines with markers represent the CSCA of the sequence obtained by the final step of the proposed algorithm. (b) Polygonal representation and areas estimation under the CSCA curves for each cluster in the ground truth shown in (a).

Fig. 5 :

 5 Fig. 5: Values of coefficients p 2 for the 1001 combinations of 4 features among 14 followed by the 2002 combinations of 5 features in binary order (points).

Fig. 6 :

 6 Fig. 6: Boxplot of NMI and IQR values against the number of clusters for all simulated datasets. The central mark is the median, the edges of the box are the 25-th and 75-th percentiles, the whiskers extend to the most extreme data points, and the outliers are plotted individually (red crosses).

  (a) Uncertainty on cluster 1 before fusion. (b) Uncertainty on cluster 1 after fusion. (c) Uncertainty on cluster 2 before fusion. (d) Uncertainty on cluster 2 after fusion. (e) Uncertainty on cluster 3 before fusion. (f) Uncertainty on cluster 3 after fusion.

  Figure8illustrates the performances over all subsets using a box plot representation for the three algorithms. On average, the best performance (96%) was obtained for the proposed method using the GK algorithm. When using the HMM, the performance in terms of precision was lower (90%) as well as for the robustness (interquartile range equals 22% for HMM against 9% for GK). The Kmeans-based method led to 12%. The important difference between the two methods (voting based on Kmeans with single feature against multiple clusterings) in terms of performance holds in the fact that the Kmeans optimises the compacity and the circularity of clusters[START_REF] Wu | Top 10 algorithms in data mining[END_REF] which is not relevant for this dataset. On the opposite, the proposed method looks for sequence of clusters. The use of the GK algorithm allows to coping with AE data scattering, where clusters can take the form of ellipsoids, through the use of a Mahalanobis-like distance[START_REF] Gustafson | Fuzzy clustering with a fuzzy covariance matrix[END_REF]. In addition, the Mahalanobis distance is known to be robust to outliers[START_REF] Peel | Robust mixture modelling using the tdistribution[END_REF]. The proposed method is thus more appropriate when the dynamics follow a sequence of states, when the distribution of features in each state is not globular-shaped and when the number of elements in each cluster is different. It is important to point out that some individual sequences obtained by some subsets of features led to erroneous sequence using the GK clustering algorithm as shown in Figure4and in Figures7(a), 7(c) and 7(e). However, those errors were compensated by the fusion process that draws benefits from multiple clusterings (Figure7(b), 7(d) and 7(f)).

  (a) IR images emphasizing hoop splitting of 90 o layers and delamination. The points 1 . . . 4 represents the location of the measurements to establish the evolution of the temperature (depicted in Fig. 11(b)). (b) Evolution of the relative increase in temperature at the four points depicted in Fig. 11(a), hoop stress and cumulated acoustic energy as a function of loading time.

Fig. 11 :

 11 Fig. 11: IR images obtained with a Mid-Wavelength InfraRed camera with 3-5 m spectral response (MWIR3 JADE from CEDIP Company).

  Cluster fusion (bold lines) and individual partitions (gray level).

  n a l s t r e n g t h A b s o l u t e e n e r g y (b) Feature representativeness in the selected individual partitions used in the fusion process.

Fig. 10 :

 10 Fig. 10: Distribution of amplitudes in clusters estimated after fusion, and feature representativeness in the fusion process.

Fig. 12 :

 12 Fig. 12: Optical observation of the edge just before the total ruine, showing different damage types (obtained by a CCD Ueye 1MPixel camera with a frame rate of 1 frame per second): A) Ply delamination, B) 45 o fibres failure, C) 90 o fibres failure.

Fig. 13 :

 13 Fig. 13: Fracture facies (using an optical microscope Nikon Eclipse LV 150) after specimen failure: 1) Hoop splitting of the 90 layers at the edge of the notch; 2) Ply delamination; 3) Fibres breakage.

ACKNOWLEDGMENT This work has been carried out in the framework of the Laboratory of Excellence ACTION through the program "Investments for the future" managed by the National Agency for Research (references ANR-11-LABX-01-01). The authors are grateful to the Région Franche-Comté and "Bpifrance financement" supporting the SMART COMPOSITES Project in the framework of FRI2. The authors wish to express their gratitude to M.G.R. Sause to share his dataset.