Progress in Aerospace Sciences A new methodology for automating acoustic emission detection of metallic fatigue fractures in highly demanding aerospace environments: An overview

The acoustic emission (AE) phenomenon has many attributes that make it desirable as a structural health monitoring or non-destructive testing technique, including the capability to continuously and globally monitor large structures using a sparse sensor array and with no dependency on defect size. However, AE monitoring is yet to ful ﬁ l its true potential, due mainly to limitations in location accuracy and signal characterisation that often arise in complex structures with high levels of background noise. Furthermore, the technique has been criticised for a lack of quantitative results and the large amount of operator interpretation required during data analysis. This paper begins by introducing the challenges faced in developing an AE based structural health monitoring system and then gives a review of previous progress made in addresing these challenges. Subsequently an overview of a novel methodology for automatic detection of fatigue fractures in complex geometries and noisy environments is presented, which combines a number of signal processing techniques to address the current limitations of AE monitoring. The technique was developed for monitoring metallic landing gear components during pre- ﬂ ight certi ﬁ cation testing and results are presented from a full-scale steel landing gear component undergoing fatigue loading. Fracture onset was successfully identify automatically at 49,000 fatigue cycles prior to ﬁ nal failure (validated by the use of dye penetrant inspection) and the fracture position was located to within 10 mm of the actual location. a variety of mechanisms, such as crack tip advance, plastic deformation, or other mechanical behaviour like friction and rubbing. This energy radiates from its point of release, known as the source, in all directions, propagating as an elastic stress wave. Minute surface displacements generated when a wave is incident upon the surface of a structure are detected using piezoelectric transducers, the voltage response of which is stored digitally for analysis. The nature of this technique means that a source mechanism must be active (i.e. damage must be growing) in order for it to be detected, making it ideal for in-service structural health monitoring (SHM) and non-destructive testing (NDT). The technique allows global monitoring and localisation of the source position within large-scale structures using a distributed array of sensors and, given that detection requires only for a source mechanism to be active, it is not dependent on defect size. Despite this potential, poor performance of the AE technique in complex geometries and high background noise levels found in industrial environments has limited application. A number of challenges related to accurate localisation and classi ﬁ cation of sources still remain, which must be addressed in order to exploit the full potential of the AE technique. This paper presents an overview of a new methodology for accurate and automated location and identi ﬁ cation of fatigue cracking in complex geometry, high noise, industrial environments using the AE technique. This represents the culmination of an extensive research program in which experimental practices and signal processing techniques have been developed to address the limitations of the AE technique. The methodology was developed for the detection of fatigue cracks during pre- ﬂ ight certi ﬁ cation testing of aircraft landing gear components; an industrial fatigue environment where complex geometries, multiple moving parts, varied load scenarios (requiring up to


Introduction
Acoustic Emission (AE) is the term given to the physical phenomenon whereby small amounts of elastic energy are released within a structure by a mechanical mechanism. Such energy release may arise from a variety of mechanisms, such as crack tip advance, plastic deformation, or other mechanical behaviour like friction and rubbing. This energy radiates from its point of release, known as the source, in all directions, propagating as an elastic stress wave. Minute surface displacements generated when a wave is incident upon the surface of a structure are detected using piezoelectric transducers, the voltage response of which is stored digitally for analysis. The nature of this technique means that a source mechanism must be active (i.e. damage must be growing) in order for it to be detected, making it ideal for inservice structural health monitoring (SHM) and non-destructive testing (NDT). The technique allows global monitoring and localisation of the source position within large-scale structures using a distributed array of sensors and, given that detection requires only for a source mechanism to be active, it is not dependent on defect size. Despite this potential, poor performance of the AE technique in complex geometries and high background noise levels found in industrial environments has limited application. A number of challenges related to accurate localisation and classification of sources still remain, which must be addressed in order to exploit the full potential of the AE technique.
This paper presents an overview of a new methodology for accurate and automated location and identification of fatigue cracking in complex geometry, high noise, industrial environments using the AE technique. This represents the culmination of an extensive research program in which experimental practices and signal processing techniques have been developed to address the limitations of the AE technique. The methodology was developed for the detection of fatigue cracks during pre-flight certification testing of aircraft landing gear components; an industrial fatigue environment where complex geometries, multiple moving parts, varied load scenarios (requiring up to 30 actuators) and very high levels of background noise are common. These test can last up to four years, with a significant proportion allocated to periodic non-destructive testing (NDT), hence the potential cost savings from implementing an online monitoring tool are large. The paper begins by outlining the challenges that must be overcome in the development of such a methodolgoy and reviews the progress made in addressing these challenges to date. An overview of the new methodology is given and previously published contributing work is referenced where appropriate. Finally the results of validation testing are presented from a full-scale landing gear component undergoing fatigue loading.

The challenges
In outlining the challenges of accurate location and characterisation of AE signals it is appropriate to first consider the AE transfer function (the relationship between source mechanism and digitised signal), presented graphically in Fig. 1, the solution of which is a complex and highly uncertain inverse problem. It shows that the form of the digitised AE signal depends upon the wave propagation path from the source to the sensor (i.e. material and geometry), the coupling arrangement of the sensor, the transfer function of individual sensors and the acquisition hardware. Thus changing any of these dependancies will result in a different form of recorded signal from an identical source mechanism.

Signal descriptors
Although analysis of AE signals can be performed using the full digitised transients, this requires a large amount of data storage and computational effort. In practice it is common to use a set of features that describe the form of a signal. Traditionally these are time based features such as signal amplitude, signal length (duration) and time to reach peak (rise time) [1], however, signal features derived in this way have several weaknesses. Since the parameters are dependent on a user defined threshold they are not universally comparable, furthermore many features are interrelated, e.g. as the amplitude of a signal increases so the energy necessarily increases. Finally, and perhaps most importantly, time based feature extraction does not capture the underlying modal structure of an AE signal [2]. Effective signal descriptors must describe the detailed form of a signal using the least possible data.

Signal classification
The AE transfer function shown in Fig. 1 highlights the complexity of the classification task using AE. In the case of a single sensor mounted to a complex structure (where couplant, sensor and system characteristics remain unchanged), the transfer function will vary with position across the structure, due to varying propagation paths. For the case of multiple sensors on a complex geometry, a single source mechanism will have a differing propagation path to each sensor, each with a differing transfer functions. Therefore an array of n sensors will acquire n different signals from a single source mechanism in a given position. This problem grows increasingly complex when variation in the mechanics and orientation of a given source type (i.e. crack growth) is considered.

Location
The capability to locate damage is one of the most powerful attributes of the AE technique and traditionally this is achieved using the time of arrival (TOA) algorithm, detailed in the NDT Handbook [3]. This is the approach adopted by current commercial systems and a location is resolved through the minimisation of the objective function in Eq. (1) by adjusting the x and y positions: where t Δ i obs , are experimentally measured wave arrival time differences between each sensor pair combination and t Δ i calc , are calculated arrival time differences, based on an assumed source position (x,y), the known sensor positions and a single user specified wave velocity. The assumed source position that gives the minimum output of Eq. (1) is taken to be the damage location. The limitations of this approach are discussed in detail in the NDT Handbook [3] and by Eaton et al. [4] and fall in to two main categories: those related to signal arrival time measurement and those related to inaccurate representation of the propagation path and velocity. The arrival time is traditionally taken as the point at which the signal amplitude exceeds a pre-set threshold level, however, multiple signal peaks often exist prior to this point. The difference between the true and measured signal arrival is one source of error and a second arises from the frequency dependent propagation velocity which requires that the signal arrival be triggered by the same phase at each sensor to ensure accuracy. In the TOA approach, signal propagation is assumed to be within a 2-dimensional plane, the source to sensor path is assumed to be direct and uninterrupted and the propagation velocity is assumed to be constant in all directions and regions of the structure. In practice this is rarely the case with curved surfaces, holes, lugs and thickness changes rendering these assumption invalid and contributing to a reduction in accuracy in complex materials and geometries.

Previous progress
A significant body of research work has focussed on addressing the limitations of AE technique discussed above and the following is a review of the current state of the art.

Classification
Early attempts at classification using a single feature, such as amplitude or peak frequency [5][6][7], or 2D correlation plots of two features [1] are now generally considered insufficient to reliably separate signals from different sources. To ensure reliable classification a greater number of signal features is required and hence requires the adoption of more complex multivariate pattern recognition techniques. This type of classification falls in to one of two categories: supervised or unsupervised. Supervised classifiaction requires the use of a suitable a priori set of training data, whereas unsupervised classification aims to group signals based on their similarity with no a priori information. A range of approaches to supervised classification have been explored, which include artificial neural networks (ANN) [8][9][10] and Support Vector Machines (SVM) [11,12]. However, the dependence of signal features on the position of the source, means an extensive training set is required containing each damage mechanism from positions across the whole structure. In the case of unsupervised classificaiton a range of common clustering algorithms have been widely adopted that include k-means and fuzzy c-means [13][14][15][16] and Pomponi and Vinogradov resented a real-time implementation of the k-means algorithm [17]. Successful implementations of both supervised and unsupervised classification approaches have focussed on small coupon scale samples where the effects of propagation, geometry and sensor transfer function are negated, but of course extrapolation up to an industrial scale is then highly problematic. Rippengill et al [18] used a number of automatic classification approaches, including Gaussian statistical methods, kernel density estimation and multi-layer perceptron ANNs, to separate noise and fatigue cracking signals from a steel bridge box girder. Crucially, in this large-scale structure, they discovered that data from near field fatigue cracking and far field fatigue cracking were also separated, highlighting the effects of propagation on classification in large structures. Al-Jumaili et al [19] also demonstrated very poor classification of artificial signals in large composite structures and went on to develop the parameter correction technique (PCT) which creates a set of propagation corrected signal features that allowed accurate grouping of similar artificially generated signals. The task of attributing clusters to damage mechanisms still remains a challenge.
As an alternative to pattern recognition, a number of researchers have exploited the modal nature of AE waves in plate like structures, using the amplitudes of the s 0 and a 0 modes to determine source orientation in both composite materials and steels structures [20][21][22]. Other researchers have developed approaches to account for the effects of propagation on the mode amplitudes [23,24] and therefore improving accuracy in larger structures. However, in all cases the approach relies on the successful identification of the two primary propagation modes, which is a challenging process to automate and the presence of reflections from boundaries makes this challenge far greater [25].
This dependence on the effects of propagation has inspired great interest in the prediction of wave propagation behaviour. Standard numerical solutions to Lamb's wave equations, although informative, are not sufficient for this task as they do not account for attenuation and damping within the material, nor do they account for the effects of boundaries and structural features. What is needed is a full-field 3dimensional model of wave propagation. This has been achieved using a range of approaches including local interaction simulation approach 'LISA' [26,27], Finite Element analysis [28][29][30][31] and spectral element analysis [32]. Despite achieving some accurate results, such modelling approaches require significant computational power and are not practical for classification in an industrial environment where signal numbers can reach 100,000 s per day. Wilcox et al [33,34] proposed a modular 'forward model' transfer function approach in which each stage of the AE transfer function, is modelled and has been expanded to include anisotropic materials [35]. The approach performs well in simple parallel sided plates, but it is suggested that further increases in complexity will require the implementation of an FE based propagation model.

Location
As outlined in the challenges above, reliable wave arrival time measurement is paramount to achieving accurate calculation of AE source location. The traditional threshold crossing technique, although simple, is ineffective and doesn't give repeatible results. A range of approaches to determine the arrival of specific frequencies or propagating modes have been investigated. These include narrow band filtering and time frequency decomposition techniques such as wavelet transformations [1,[36][37][38][39][40], however, these techniques become less accurate in the presence of wave reflections from boundaries and structural features [37]. A prefered approach is to detect the first signal motion above background noise. The use of statistical approaches to separate random noise prior to signal onset and structured signal following signal onset have shown much promise in pursuit of this.
Lokajicek and Klima [41] take the derivative of the 6th order statistical moment of data within a sliding window; the normal distribution of noise data is skewed with the inclusion of a few signal points in the window, thus affecting a change in the 6th order moment. Up to 95% of signal arrivals were estimated to within 2 samples (at 10 MHz sampling rate) using this method. The use of a threshold is still required to detect the corresponding change, which can be a further source of error. Others [42][43][44] have utilised an approach based on the Akaike Information Criterion (AIC) [45] to determine signal arrival times. The AIC function compares signal entropy before and after each time t (or data point) in a time series. When time t is aligned with the signal arrival, the similarity between the high entropy uncorrelated noise prior to t and the low entropy structured signal after t is at its lowest and the function returns a minimum. A simple minimum finding function can then be used to reliably determine the signal arrival. This approach has been shown to be very reliable with arrival time estimations varying between 2-4% when compared with manually selected arrival times. The AIC approach is adopted in the methodology presented below and will be discussed in greater detail.
To account for the anisotropic propagation velocity commonly observed in composite materials Paget et al [46] developed a closed form solution based on the assumption of an elliptical wave front. The assumption of an eliptical wave front is not valid in many composites and closed form solutions are rarely stable enough when considering real test data. A more appropriate solution updates the TOA algorithm with a velocity term that references a user defined velocity profile representing the material in question [47,48]. As an alternative Ciampa and Meo [49] proposed a novel triangular arrangement of closely spaced sensor pairs (six sensors in total) for which six non-linear equations with 6 unknowns (source position, x and y, travel time to first sensor and three propagation velocities) can be derived. Solving iteratively using a Newton method provides a source position without the need for a priori knowledge of the velocity profile. Although effective, the algorithm is computationally demanding requiring 2 s to resolve each location and required sensor numbers are doubled for a given area. Despite these promising results none of these approaches are capable of addressing the structural complexity experienced in an industrial environment. All demonstrations have been in simple geometries (flat plates) where propagation paths are linear and velocity does not change along the path. In reality propagation paths will interact with structural features (holes, lugs etc.) and velocity will change regionally as thickness changes or in composites as the layup changes. Horn [50] went some way to addressing this with a reverse ray tracing approach in which reflections from boundaries etc. are considered. However the model is computationally expensive requiring large amounts of parallel processing, hence not practical for real time monitoring applications.

Methodology
This section outlines the proposed methodology of this work, in which accurate detection and location of fatigue fracture can be achieved in highly complex and noisy structures, automatically and without the need for operator input. The process flow of the methodology is outlined in Fig. 2 and each step is discussed in more detail below.

Arrival time and feature extraction
The importance of arrival time estimation is discussed above and in this methodology an AIC [45] based approach has been adopted in order to provide accurate, threshold-independent arrival time estimation. To achieve this signal acquisition is triggered using the traditional threshold approach and arrival times are subsequently corrected using the AIC appraoch, as discussed by Hensman et al [42]. For the purposes of this paper a short description of the AIC index is given: 10 10 where var is the classic variance of a given vector. As can be seen, the AIC functions in a direct way by splitting the signal into two vectors one given by x t { (1: )} for time 1 to t and a second vector x t T { ( : )} from time t to T. The core of the AIC (or any other information-based criterion) is to characterise the entropy (i.e.information context) by describing the similarity between the two broken parts of the total signal in such a way that the part before t, includes high-entropy context with uncorrelated noise and the part after t in contrast contains low-entropy context of given correlation. The target is to estimate the minimum of the AIC at the later point. Fig. 3 demonstrates the use of the first threshold crossing arrival time estimation (grey lines) and the AIC arrival time estimation (black lines), where a clear improvement is seen.
To ensure accurate classification, signal descriptors (or features) must be selected that most accurately represent the form and the true modal nature of an AE signal and this can be achieved by considering signals in the time-frequency domain. Initially, the use of Gaussian Mixtures Modelling was adopted to describe the distribution of energy within a time-frequency wavelet transform of the signal. This was shown to perform favourably when compared with traditional signal features for the classification of AE data [51]. Subsequently a Fast Wavelet Transform (FWT) algorithm was developed [51] that provides much greater computational efficiency when decomposing a signal in time-frequency, whilst still providing an accurate representation of the signal form. The decomposition process is repeated until a small number of vectors remain, 32 in the case of this work, which describe the original signal. These signal features can then be stored at very little storage cost for subsequent processing.
Generally, a wavelet transform gives the time-scale response of a signal with the result being a surface, with the axes representing time, scale or frequency, and amplitude of the analysed signal. The wavelet decomposition of the signal in this work is based on deriving subspaces of the signals by approximating the wavelet multiple scales by a series of Gaussians: where A n is the height of the nth Gaussian, x o n , is the horizontal position of the nth Gaussian, y o n , is the vertical position of the nth Gaussian, σ x n , is the horizontal width of the nth Gaussian and σ y n , is the vertical width of the nth Gaussian.
Each Gaussian presents a decomposed vector in the surface. The selection of a suitable number of decomposed vectors is a key part of the process as with only a few of them, the approximation of the wavelet transform is poor, and the subspace will not accurately represent the total signal. On the other hand with too many of them, the subspace will become very large, and not practical as a tool.
There are multiple methods that one could use to optimise the position and sizes of each Gaussian and the key point is to automatically determine the position and size each Gaussian so it can accurately and efficiently approximate the wavelet surface of the whole signal. Methods such as MCMC or evolutionary algorithms (like genetic algorithms or differential evolution) or fast simulated annealing can be adopted. For the purposes of the current work the reader is referred to [51].

Source location
Source location forms the 'front end' of the proposed methodology and as such it is essential to ensure accuracy and reliability, even in complex material and geometries. To achieve this Baxter et al [52], developed an empirical approach to account for high levels of structural complexity, known as Delta T Mapping. The same approach has shown excellent performance in composite materials [53] and is outlined in the following steps: • Construct a map system on area of interest A grid is placed over the area of interest within which AE events will be located and locations will be calculated with reference to the grid coordinates.
• Obtain time of arrival data from an artificial source An artificial source (such as a Hsu-Neilsen source [54,55]) is generated at grid nodes to provide AE data for each sensor. An average of several repeated sources is used for each node and missing data points can be interpolated.
• Calculate the DeltaT map Each artificial source results in a difference in arrival time or Delta T for each sensor pair (an array of four sensors has six sensor pairs). The average Delta T at each node is stored in a map for each sensor pair. The resulting maps can be visualised as contours of constant DeltaT.
• Locating real AE data The DeltaT values from a real AE event are calculated for each sensor pair. A line of constant DeltaT equivalent to that of the real AE event can then be identified on the DeltaT map and by overlaying the resulting contours for each sensor pair, a convergence point is found at the source position. In theory, all the lines should intersect at one location, however in practice this is not the case. Thus a cluster analysis is used to estimate the most likely source position.  4.2.1. Gaussian process regression algorithm Hensman et al. [42] extended this concept through the use of Gaussian Processes (GPs) to learn the relationship between the arrival times at a group of sensors and the corresponding positions on the structure. The GP regression was used to map the vector T Δ for each pair of sensors positioned on the experimental structure. The AIC arrival time estimation algorithm was also included in this approach. This work offered a number of improvements over the original DeltaT Mapping algorithm: reduced training data, increased computational efficiency, a measure of confidence due to the probabilistic approach, the ability to generalise across similar structures and the automatic selection of the best sensors. AIC was very important in terms of using a Bayesian regression system like GPs as the signal onset can influence the performance of GPs when used to locates AE sources. A short description of GPs is given below as it would be helpful for the reader.
Rasmussen and Williams [56] define a Gaussian process (GP) as "a collection of random variables, any finite number of which have a joint Gaussian distribution". In recent years, GPs are gaining a lot of attention as a machine learning approach in the area of regression (or classification) analysis as they offer fast and simple computations. Gaussian process regression is a robust tool which takes into account all possible functions that fit to the training data set and gives a predictive distribution rather than a single prediction for a given input vector. As a result, a mean prediction and confidence intervals can be calculated from this predictive distribution. The basic details of the algorithm are presented following the steps in [56]. The algorithm that was used in the previous sections is also coming from Rasmussen and Williams [56].

Algorithm theory
The initial step in order to apply Gaussian process regression is to define a prior mean m where E represents the expectation. Often, for practical reasons because of notation purposes (simplicity) and little knowledge about the data at the initial stage the prior mean function is set to zero. The Gaussian processes can then be defined as: Assuming a zero-mean function, the covariance function could be described as: This is the squared-exponential covariance function (although not the only option). It is very important to mention an advantage of the previous equation as the covariance is written as a function only of the inputs. For the squared-exponential covariance, it can be noted that it takes unit values between variables where their inputs are very close and starts to decrease as the variable distance in the input space increases.
Assuming now that one has a set of training outputs f { } and a set of test outputs f { } * one has the prior: where the capital letters represent matrices. As can be seen, the covariance matrix must be symmetrical about the main diagonal.
As the prior has been generated by the mean and covariance functions, in order to specify the posterior distribution over the functions, one needs to limit the prior distribution in such a way that includes only these functions that agree with actual data points. An obvious way to do that is by generating functions from the prior and selecting only the ones that agree with the actual points. Of course, this is not a realistic way of doing it as it would consume a lot of computational power. In a probabilistic manner, the operation can be done easily via conditioning the joint prior on the observations and this will give (for more details see [56][57][58]): Function values f { } * can be generated by sampling from the joint posterior distribution and at the same time evaluating the mean and covariance matrices from (9).
The covariance functions used in this study are usually controlled by some parameters in order to obtain a better control over the types of functions that are considered for the inference. As an example, a more general form of the squared-exponential covariance function can take the form (1-dimensional): where k y is the covariance for the noisy target set y (i.e.y f x ε = ({ }) + where x { } is input vector and ε is the noise). The length-scale l (determines how far one needs to move in input space for the function values to become uncorrelated), the variance σ f 2 of the signal and the noise variance σ n 2 are free parameters that can be varied. These free parameters are called hyperparameters.
The tool that has to be applied for choosing the optimal hyperparameters for GP regression, is the maximum marginal likelihood of the predictions θ p({y}|[X], { }) with respect to the hyperparameters θ: [ ] f is the noise-free covariance matrix. In order to optimise these hyperparameters through maximising the marginal log likelihood the partial derivatives give the solution, via gradient descent: Of course this solution is not a trivial procedure and for specific details readers are refereed to [56].

Spatial clustering and fracture detection
Spatial clustering groups together signals origniating from the same location, and hence the same source mechanism. It allows properties of a cluster to be considered as opposed to individual signals, giving greater statistical relevance. This is an important concept for this methodology and further highlights the importance of accurate source location. Spatial clustering has been implemented through the segmentation of the area of interest using a regular grid and also through the use of an on-line radius-based clustering algorithm (ORACAL) [59]. The ORACAL approach is a data driven clustering approach that allows the merging and separation of clusters as the data set evolves, making it ideal for a real time application such as this.
Fracture detection is achieved through the consideration of two cluster properties: the normalised signal feature variance and the novelty of spatial location. These properties are chosen based on the following two hypotheses: • A fatigue crack is a repeatable mechanism that will produce repeatable signals, whereas noise sources, e.g. friction and rubbing, are more random and produce a greater variation in signals.
• Noise sources will always be active and hence produce a distinctive location pattern or "footprint". Fracture onset results in a new source mechanism, causing a change in the established footprint.
Variance of signal features within a cluster is considered on a perchannel basis and an example of this is presented in Fig. 4. All shaded clusters contain located AE and all three sensors are required to locate these data, however to assess the variance feature within a cluster only one of the three sensor is used. In this example the variance of the red cluster is considered using channel 1, hence all the data have the same propagation path, coupling arrangement and sensor transfer function. Therefore any variation observed can be attributed to the variance in the source mechanism. The variance of a cluster is also normalised against all other data located by the considered channel. In this way, the lower variance of a cluster resulting from a fatigue fracture should be further highlighted with respect to higher variance noise data. The process can also be conducted for all other sensors but only data from one sensor can be considered at any one time and comparison between different sensors is not considered valid. This condition is an important feature of the analysis methodology and negates the need for computationally intensive classification. The use of cluster variance has been successfully demonstrated for the detection and identification of artificially generated fatigue fracture signals in an industrial fatigue environment [59]. The computation of multi-dimensional variance requires that a cluster contains an equal or greater number of signals than there are dimensions, or number of signal features used, which in this case is 32. In the analysis below a cluster population is required to exceed 40 signals before a variance value is computed. The resulting normalised variance is then presented on an inverse logarithmic scale.
The identification of the spatial novelty of a cluster location utilises the spatial scanning technique, proposed by Hensman et al [60], that identifies the appearance of a new source against an established background. A summary of the proposed method is described below: First of all, it has to be clear that spatial scan statistics are performing an extensive scan of the spatial distribution of the data in order to identify areas with noticeable activity, or in other words it is a clustering method [61,62].
If one assumes a series of locations within the data cloud that are indexed by i then each location has a population b i which could be at risk q of presenting novelty (or a disease symptom as the method is coming from 'syndromic surveillance'). If a fixed time window t of the given data is assumed then the appearance of novelty at any location c i could follow a Poisson distribution: The Poisson distribution PO is suitable for AE emission data as in this case one talks about big data where each population within this data is considered to be at a small risk. In order to practically perform scanning of the online data for over-activity, one must introduce a definition of all potential clusters. If S is the set of all indices on different locations and K is a subset of the indices S (which can be defined as S k for k K = 1, 2, …, ) which represents the data space one wishes to inspect as a cluster then the joint-likelihood of all data C given by K sets of clusters with given indices can be written as: where S k c is the rest of the space S k , B is the sum of all baselines b i and q in ,q out are the risk of novelty within a cluster and outside of the specific cluster respectively. In the AE methodology the baselines b i can be calculated as: which represents the expected number counts at position i for baselines b i at time t which represents the global activity duration. The remaining thing is the calculation of priors q in and q out which can be implemented by considering previous data that does not exist within a cluster. In such case S k =0 and the likelihood of data p C B α β ( | , , ) can be computationally simplified. A Γ distribution can generally be defined as: where c is any positive real number.
x is a multiplier when there is a risk increase by the indication of a cluster, α β , are the prior parameters of a Γ distribution that can be obtained by utilising a Markov chain Monte Carlo (MCMC) tool. In order to get a valid hypothesis test for values q in and q out then the last equation has to be marginalised. Lastly, x can be obtained by prespecified values given by [61,62] with equal probabilities each of them.
A statistical representation of the spatial distribution and activity rate of recent test data is compared with that of previous data from the test. Any areas of abnormally high activity observed in the recent data will register as an anomaly within the distribution, indicating the occurrence of a new source, such as the onset of a fatigue fracture. The spatial scan statistic does not attempt to perform source characterisation: it bypasses this concept by simply looking for abnormal levels of activity. It is a novelty detector. Cluster activity rate has also been considered as a fracture indicator, under the assumption that a noise source should give an approximately linear activity rate whereas a  fracture source would give an increasing activity rate in line with the Paris-Erdogan law for fatigue. However this was found to be an unreliable indicator and as such, is no longer used for fracture identification.

Visualisation and alarm
Appropriate visualisation is an important step in the proposed methodology to ensure the best representation of the analysed data. Data are initially viewed as location plots, where clusters containing AE events are shaded grey (a darker shade means more events) and those containing greater than 40 events are assigned a coloured dot, the colour of which indicates the cluster variance score (see Fig. 8). The variance score is given on an inverse log scale, hence larger numbers indicate fracture. Plots such as this can become saturated with data in long duration, high noise tests, such as landing gear certification testing. So for the purposes of clarity only data occurring within the preceding t seconds are considered in this view, however cluster variance is normalised against all data recorded up to this point. Similarly, the novelty data can be viewed as a location plot where clusters are shaded if they vary from the established "footprint" (see Fig. 10). Again, data is taken from the preceding t seconds of testing and compared with the "footprint".
For the purposes of alarm, trends of maximum observed spatial novelty and maximum observed cluster variance across all sensors are considered. The earlier stages of testing can be used to provide a baseline level of values against which a threshold could be set to raise an alarm when values increase. In order to reduce the likelihood of any false positives, fracture identification is based on a simultaneous activation of both indicators. Following alarm an operator can perform a more detailed analysis identifying the location of the clusters. If the same cluster location is identified by spatial novelty and signal variance then the user can be satisfied that a confident identification and location of a fracture has occurred.

Validation testing
In order to validate the presented methodology, a steel A320 landing gear main fitting was subjected to fatigue loading. Fig. 5 presents an image of the component and shows the high geometric complexity. A fatigue crack source was generated by cyclically loading a representative steel lug welded at the site of a previously tested lug. Although welding reduces the strength of the lug compared with the original component, the validation testing requires only that a fatigue fracture be generated. The lug, shown in Fig. 6a, is 40 x 25 mm and 10 mm thick and was loaded under bending using the cantilever loading arm shown in Fig. 6b. The lug was loaded at a frequency of 1 Hz and with an R ratio of 0.1, initially with a peak load (applied vertically to the loading arm) of 5.5 kN, then at 6 kN after 90k cycles, 6.5 kN after 110k cycles and 7 kN after 138.5k cycles until fatigue fracture at 168.7k cycles. To replicate the high levels of noise experienced during landing gear certification testing the sliding tube, seen in Fig. 5 was also activated throughout testing at a frequency of 0.2 Hz. Loading was stopped periodically (every 5k cycles) throughout the test, the loading arm removed and a dye penetrant inspection performed on the lug and weld to identify the onset of fracture. The main fitting was monitored continuously using eight Mistras Group Limited (MGL) Nano30 sensors, arranged in an array around the cylindrical part of the landing gear component. The sensor positions and main geometric features of the component are shown schematically in Fig. 7, for the area of interest. The sensors were attached to the component using magnetic clamps and acoustically coupled using brown grease. The sensitivity of the installed sensors was determined using a H-N source, in accordance with ASTM standard E976 [54,55] and the response of all installed sensors was greater than 97 dB. All AE data were recorded using a Physical Acoustics Ltd. PCI-2 acquisition system and waveforms for all signals were sampled at 2 MHz throughout the investigation. The collection of waveforms was triggered using threshold crossing (43 dB) and a hit lock out time of 250 μs (which specifies the time that a sensor output should remain below the threshold before any subsequent threshold crossing will be considered as a new hit).
Data analysis was conducted post-test using the stored sampled waveforms. An area of interest of 0.75 x 0.7 m was divided into spatial bins of 19 x 18 mm to form clusters for both novelty detection and variance calculation. In application to real-time monitoring the data of most interest is that which occurs within the last "t" seconds and allowing data to "age" is important to ensure detection is not compromised. To replicate real-time analysis in this way in a post processing scenario, a sliding time window of duration t seconds is used. Selecting the length of time window within which to interrogate the data is very important and it may be that an optimum time window can be determined from the expected data rates and monitoring duration. However, for the purposes of demonstrating this methodology a trial and error approach is sufficient and time windows of t=1800 s and t=3000 s were selected for novelty detection and variance score calculation, respectively.

Results
Dye penetrant inspection of the loaded lug revealed the first possible indication of fracture, with low confidence, after 125,000 s of testing and the indication remained unchanged after 130,000 s (at 1 Hz loading frequency seconds are equivalent to cycle number). Confident detection of the fracture was observed during the dye penetrant inspection after 135,000 s of testing and final failure occurred after 169,000 s Fig. 8 presents the located clusters within the area of interest after 60,000 s of testing, 65,000 cycles prior to the first indication of fracture from dye penetrant. The key features of the plot will be discussed prior to consideration of the results. The boxed numbers represent the position of sensors and the red coloured sensor number indicates the "active sensor"; therefore only events located using this sensor are presented and only signals recorded by this sensor are used for the calculation of the presented cluster variance scores. The grey rectangles indicate spatial bins (or clusters) that contain 1 or more events. Clusters containing more than 40 events are attributed a variance score indicated by a coloured circle. A colour key on the right hand side of the plot represents the value of calculated variance scores. It is important to note that the data presented in the location plot are from the previous 3,000 s of testing (i.e.the time window, 57,000-60,000 s), however the variance score for a cluster is normalised against the entire data set (i.e. all located signals recorded by the active channel for 0-60,000 s) up to that point. The higher the  normalised variance score the higher the likelihood that the cluster is from a fracture source, acknowledging that the variance scores are presented on an inverse logarithmic scale. Finally, the red cross in the location plots indicate the central position of the tested lug at 0.33, 0.05. A number of clusters containing more than 40 events are identified and attributed variance scores of up to 70 and little activity is observed in the region of the lug. It is not possible at this stage to say if these are high or low variance scores, however, being so far prior to fracture identification using dye penetrant it can be assumed that all clusters originate from noise sources. Given data from a greater number of tests it would be possible to estimate a baseline level of variance scores, above which a cluster can be considered as originating from a fracture. None of the located clusters at this stage of the test were observed to have any spatial novelty. Fig. 9 shows a similar plot after 122,000 s of testing, 3000 s prior to the first indication of fracture onset by dye penetrant. At this stage of the test, clusters have been located directly adjacent to the lug position, one of which has exceeded the cluster threshold and has been attributed a variance score of 120. This is much higher than the values attributed to any other cluster located elsewhere in the area of interest and of those previously seen in Fig. 8, indicating that fracture growth is occurring within the loaded lug. The cluster size used in this analysis is 19 mm×18 mm, making the centre of the identified cluster ∼8 mm from the corner of the lug at which dye penetrant detection of the fracture was subsequently observed. Similar results were observed when using channels 3, 5 and 8 as active channels, these being the other channels used in the location of fracture signals. The spatial novelty of the located clusters at this stage of the test is presented in Fig. 10, where each cluster is shaded based on its level of novelty. The centre of the loaded lug, at 0.33, 0.05 m, is indicated by the grey cross and there are four clusters exhibiting spatial novelty located adjacent to the lug position and their spatial novelty is high. The combination of the high variance score and high novelty attributed to the clusters located adjacent to the lug position gives a very confident indication of fracture onset. This is further supported by the first indication of fracture by the subsequent dye penetrant inspection at 125,000 s. It is noted that detection by dye penetrant is only with low confidence at this stage, unlike the AE results.
To simplify the detection of fracture onset and to limit the operator input required during monitoring of a long-term test, the results can be viewed as trends of novelty level and variance score with time of test. Fig. 11 presents an example of this for the current test, where trends of maximum observed spatial novelty (top) and variance score (bottom) for all located clusters across all channels can be seen. Along the top    After 120,000 s, both fracture identifiers rise sharply, confidently indentifying the onset of fracture growth 49,000 cycles prior to failure. This is supported by the first possible indication of fracture during the subsequent dye penetrant inspection at 125,000 s. Using these trends it would be easy for an operator to set a threshold for each identifier, triggering an alarm when they are exceeded. Following the dye penetrant inspection at 125,000 s, both identifiers return to normal levels and they do not simultaneously increase again until 150,000 s remaining high after that, until final failure. It is postulated that the crack growth was interrupted by the minor reconfiguration of the loading arrangement, which was necessarily removed for dye penetrant inspection. In addition to this the normalised variance score rises independently at 132,000 s following which fracture indication during the dye penetrant inspection at 135,000 s was seen to improve. Consideration of the located clusters at this point, revealed a cluster adjacent to the lug position with a high variance score, indicating that some fracture growth had occurred, however the number of signals located within the cluster was seen to be low across the 3,000 second time window. The variance score again returned to normal levels following the dye penetrant inspection. With both the variance score and the dye penetrant results indicating additional fracture growth at 132,000 s, it is of some interest that detection was not achieved by the novelty indicator. An explanation can be found through consideration of the novelty detection process. Novelty is assessed by comparing the last 1,500 s of activity in a cluster with that of an established baseline. Both the position of the cluster and its activity rate contribute to any novelty. The baseline is established from all previous test data recorded prior to the final 1,500 s under consideration. Therefore, at 132,000 s, when the variance clue is seen to increase, the novelty detection baseline includes data from fracture growth observed at 120,000 s. It is likely therefore that the novelty of the cluster resulting from fracture growth at 132,000 s was limited by the inclusion of fracture activity in the baseline data set. Indeed consideration of the novelty and variance trends for a nominally identical test (Fig. 12), but without any load interruptions for dye penetrant inspection reveals a simultaneous rise in both fracture identifiers at 40,000 s, clearly identifying fracture onset and remaining high until failure at 44,000 s.

Conclusions
A methodology for the automated detection and location of fatigue   fractures in complex geometries and high noise environments has been successfully demonstrated, overcoming a number of significant and previously unaddressed challenges. A landing gear component provided a suitably complex test environment for validation and fracture onset was clearly and automatically identified 39,000 cycles prior to final failure and located within 10 mm of the fracture position. Analysis was conducted post test in a pseudo-real-time fashion and showed that fracture onset would have been detected over 13 h prior to final failure and confident detection was achieved ahead of dye-penetrant inspection. It has been shown that through the careful application of signal processing, some of the key limitations facing the AE technique can be addressed, without the need for complex or demanding approaches such as modelling. It must be noted however that this application has focussed specifically on the automated detection of fatigue fractures in metallic components and its direct application to other materials may be limited. However, many elements of the developed methodology can offer immediate benefits, for example the demonstrated improvements in location accuracy. Furthermore the concept of spatial clustering and analysis of grouped data on a sensor by sensor basis will be a useful approach in many applications, although it is likely that new identifiers may be required.