The damage detection problem in mechanical systems, using vibration measure- ments, is commonly called Structural Health Monitoring (SHM). Many tools are able to detect damages by changes in the vibration pattern, mainly, when dam- ages induce nonlinear behavior. However, a more dicult problem is to detect structural variation associated with damage, when the mechanical system has nonlinear behavior even in the reference condition. In these cases, more so- phisticated methods are required to detect if the changes in the response are based on some structural variation or changes in the vibration regime, because both can generate nonlinearities. Among the many ways to solve this problem, the use of the Volterra series has several favorable points, because they are a generalization of the linear convolution, allowing the separation of linear and nonlinear contributions by input ltering through the Volterra kernels. On the other hand, the presence of uncertainties in mechanical systems, due to noise, ge- ometric imperfections, manufacturing irregularities, environmental conditions, and others, can also change the responses, becoming more dicult the damage detection procedure. An approach based on a stochastic version of Volterra series is proposed to be used in the detection of a breathing crack in a beam vibrating in a nonlinear regime of motion, even in reference condition (without crack). The system uncertainties are simulated by the variation imposed in the linear sti ness and damping coecient. The results show, that the nonlinear analysis done, considering the high order Volterra kernels, allows the approach to detect the crack with a small propagation and probability con dence, even in the presence of uncertainties. Luis Gustavo Giacon Villani Email addresses: luis.villani@unesp.br (Luis G. G. Villani), samuel.silva13@unesp.br ( Samuel da Silva), americo@ime.uerj.br (Americo Cunha Jr.) Preprint submitted to Mechanical Systems and Signal Processing August 25, 2018 Keywords: Damage detection, uncertainties quantification, stochastic Volterra series, nonlinear dynamics. 1. Introduction The process of damage detection in mechanical, aerospace and civil sys- tems and structures, based on damage features and statistical analysis of them, is commonly called as Structural Health Monitoring (SHM) [1]. The study and application of SHM techniques are motivated by the potential life-safety and eco- nomic impact of their implementation, which contributes to the development of new approaches [2]. As it is widely known, the damage-sensitive features can be sensitive to confounding e ects such as environmental or operating conditions, temperature and humidity changes, sensor bonding conditions, and others [3]. In addition, real systems and structures are subject to uncertainties, mostly related to noise in the measurements (experimental data variation), geometric imperfections, manufacturing irregularities, environmental conditions, or even, lack of knowledge about the system physics [4, 5, 6]. These uncertainties can complicate the damage detection problem and have to be considered in the SHM applications, reducing the occurrence of false alarms [7]. To overcome this problem, the damage detection can be done based on statis- tical procedures, taking into account con dence limits to the dynamic behavior of the systems and structures. Thus, the structural monitoring can be done, usually, based on: (i) novelty detection, i. e., detection of discrepant behaviors of the structure compared with the normal condition which it is known, in most cases with the establishment of statistical thresholds to the system behavior in normal condition; (ii) classi cation, in this case all states are known and the condition of the structure or system is classi ed into groups of di erent damage stages or healthy condition; (iii) regression algorithms, in which some variable or process is monitored with con dence limits. The choice of the more interest- ing approach to be used it depends on the speci c problem and on the level of knowledge about the damage (detection, location, assessment, prediction) de- sired [8]. In this sense, a lot of methods can be used, such as linear and nonlinear Principal Component Analysis (PCA - NPCA) [9] , Extreme Value Statistics (EVS) [10], Peaks Over Threshold (POT) [11], machine learning algorithms [12], neural network [13], Bayesian approaches [14], Mahalanobis distance [15, 16], and others. In all these applications, the nal objective of the methods used are the same, di erentiate the uncertainties and variabilities of the damages. On the other hand, many existing SHM features use the nonlinear character- istics of damages to detect them. In this context, cracks [17], delamination [18] or rubbing in rotor systems [19], induce linear systems to exhibit nonlinear phe- nomena and the damage detection problem becomes a problem of nonlinearities detection. So, in this situation, the damage can be detected using harmonic distortions, coherence functions, probability density functions, correlation tests, Hilbert transform and others [20]. Thus, the crack detection in beams can be performed based on the nonlinear behavior induced by this type of damage. 2 Andreaeus et al. [21] have used the nonlinear aspects of the dynamic response of a beam under harmonic excitation to relate the nonlinear resonances to the presence and size of the crack. R ebillat et al. [22] have used cascade of Hammer- stein models to detect nonlinear damage through the monitoring of the energy of the nonlinear components in the system response. The approach was applied in a simulated breathing crack problem detection and in the monitoring of a composite plate. Lim et al. [17] have shown a reference-free fatigue crack de- tection technique, based on nonlinear ultrasonic modulation. In this last work, a sequential outlier analysis was used to identify the crack presence without re- ferring to any baseline data obtained from the intact condition of the structure. All these approaches are powerful tools used to detect the presence of fatigue crack or any other damage that induces nonlinear e ects in the system behavior. However, the SHM problem increases when the intrinsic nonlinear behavior of the systems and structures is considered because the nonlinear phenomena can be confused with damage when classical SHM techniques, based on linear metrics, are used [23]. Furthermore, many mechanical systems and structures can operate with strong nonlinear behavior, that makes them exhibit complex responses containing subharmonic and superharmonic resonances, jumps, modal interactions, bifurcation, quasi-periodicity, and possible chaos [24, 25]. In this situation, the nonlinear behavior of the system in the reference condition has, necessarily, to be considered. In this regard, many approaches, in time and frequency domain, can be used to describe the nonlinear behavior, such as Hilbert transform, Narmax Models, High-Order FRFs [26, 27], restoring force surface methods (RFS) [28, 29], har- monic balance method [30] and others. Unfortunately, these nonlinearities have many individualities, in a way that it is dicult to obtain a general model that describes all the structures of interest. Among the di erent methods for nonlinearities identi cation, the Volterra series stands out, because it is a gener- alization of the linear convolution concept, allowing the separation of the system response in linear and nonlinear components [31]. The main procedure to esti- mate the Volterra kernels is the Harmonic Probing method, that was extensively used in system identi cation problems [32, 33] and damage characterization [34]. The limitation of the approach is the dependence of a parametric model of the system. Another possible formulation is the direct use of input/output signals to estimate the Volterra kernels. Tang et al. (2010) [35] have used the classical version of the Volterra series approach to detect structural variations in a rotor, using the vibration signals of di erent points as input/output signals. Tang et al. (2010) have used the classical version of the Volterra series approach to detect structural variations in a rotor, using the vibration signals of di erent points as input/output signals. The Volterra kernels coecients and the Principal Component Analysis (PCA) were used as damages features. In this study, the use of the rst 3 kernels has shown better performance, with the capability of separate all the conditions considered. However, as mentioned by the authors, the high number of parameters to be identi ed has diculted the approach implementation using an optimization process. Similarly, Xia et al. (2016) [36] 3 have used a key kernels-PSO (KK-PSO) method to identify the Volterra model, combined with neural networks, also to detect structural variation in a rotor. The implementation has been facilitated by the model reduction obtained using the KK-PSO method. The results have shown the capability to the method to detect di erent types of structural variation. In both works, the Volterra series coecients have been used as damage sensitivity feature, but they have not considered the presence of uncertainties in the damage detection process to study the capability of the approaches to detect the damages with a probabilistic con dence. Additionally, Shiki et al. (2017) [37] have used this Volterra series formu- lation, based on input/output signals, to detect structural variations related with damage, in a beam with nonlinear behavior, even in the reference condi- tion. They have shown the use of Volterra series expanded through the Kautz functions, as a reference model to predict the nonlinear system response and to monitor the structure of interest. A statistical study of the damage index proposed was done with satisfactory results, but only the variation related to measurement noise was taken into account. The authors have not considered the possible system response variation related to the uncertainties, like environ- mental and operating conditions, and others as mentioned before. Thus, in this work, it is proposed to use a stochastic version of the Volterra series approach, expanded using random Kautz functions, as a stochastic model for damages detection based on novelty detection. The Volterra kernels are estimated in a probabilistic framework, considering the data variation, that is obtained through the variabilities in fundamental frequency and damping of the studied system. To allow the description of the response variation, the Kautz parameters are treated as random variables and the Kautz functions, used to reduce the number of terms of the Volterra kernels, as random processes. This methodology gives to the proposed model the ability to describe the system nonlinear behavior and the data variation simultaneous. To detect the presence of damage, the stochastic model obtained through the stochastic version of Volterra series is used as a set of reference models, so new models are estimated in unknown condition and compared with this set of models, through a distance-based method. The capability of the Volterra series model to separate linear and nonlinear contributions in system response is used. Two methodologies are compared, considering the linear and nonlinear components of the Volterra series. The rst methodology considers only the Volterra kernels coecients and the second considers the kernels contribution. The approaches are applied in a simulated uncertain nonlinear system that de- scribes the real behavior of a nonlinear beam, and, the damage is simulated by the propagation of a breathing crack, described as a bilinear sti ness oscil- lator. The results obtained show that the methodologies proposed are able to di erentiate the nonlinear behavior and data variation of the damage presence, with better performance than the simple linear analysis. It is showed that the use of Volterra series in this situation represents a clear advantage in damage detection process, with the improvement in the results in comparison with the linear approach. In addition, the stochastic version proposed allows the dam- 4 age detection through the approach, even in the presence of uncertainties, with probability con dence. The present paper is organized as follows. Section 2 reviews the fundamental aspects of system identi cation based on Volterra series, showing the expansion of the theory to the proposed stochastic Volterra series model, and describes the methodology used to detect damage in nonlinear systems, considering uncertain- ties of di erent natures. Section 3 presents the model of the nonlinear uncertain beam simulated and the breathing crack model used to simulate the presence and propagation of damage. Section 4 shows the results obtained through the application of the approach in the studied system. Finally, section 5 presents the main conclusions and some paths for future work. 2. The stochastic Volterra series approach for damage detection This section describes the proposed methodology for modeling nonlinear sys- tems assuming data uncertainties. The approach is based on a stochastic version of the Volterra series, expanded using Kautz functions as the orthonormal ba- sis. The sections 2.1 and 2.2 summarize the deterministic version of the Volterra series expanded using Kautz functions, which is known in the literature and it was used to detect damage in initially nonlinear systems by Shiki et al. [37]. Additionally, the sections 2.3, 2.4 and 2.5 show the expansion of the classical deterministic concept to the new random version of the approach, proposed to be used to model uncertain nonlinear systems. Finally, section 2.6 presents the two new proposed damage detection methodologies, to be used considering the random version of the Volterra series. 2.1. Deterministic Volterra series Consider a discrete-time causal nonlinear system with a single output k2 Z+7!y(k) caused by a single input k2Z+7!u(k), withZ+representing the set of nonnegative integers. Through the discrete-time Volterra series, the output of this nonlinear system can be written in the form y(k) =1X =1N11X n1=0:::N1X n=0H(n1;:::;n)Y i=1u(kni); (1) where (n1;::;n)2Z +7!H(n1;:::;n) represents the -order Volterra kernel. This formalism is very versatile, in a way that it allows one to represent di erent types of nonlinear systems using the convolution concept [31]. Also, it allows one to split the system output into a sum of linear and nonlinear contributions y(k) =y1(k)|{z} linear+y2(k) +y3(k) ++y(k)|{z} nonlinear: (2) The main drawback of the approach is the diculty of convergence of the series when a large number of terms N1;:::;Nare used. In order to reduce the number 5 of terms necessary to obtain a good approximation, the Volterra kernels Hcan be expanded into an orthonormal basis, such as the Kautz functions, that is employed in this work [38, 39]. In this way, one can write H(n1;:::;n)J1X i1=1:::JX i=1B(i1;:::;i)Y j=1 ;ij(nj); (3) whereJ1;:::;Jare the number of Kautz functions used in each orthonormal projections of the Volterra kernels, ( i1;:::;i)2Z +7!B(i1;:::;i) represents the-order Volterra kernel, expanded in the orthonormal basis, and nj2Z+7! ;ij(nj) represents the ij-th Kautz lter. Thus, using the discrete Volterra series representation of Eq. (1) and the Kautz approximation (3), the system response can be approximated by the multidimensional convolution between the orthonormal kernels Band the input signal ltered by the Kautz functions, i.e., y(k)1X =1J1X i1=1:::JX i=1B(i1;:::;i)Y j=1l;ij(k); (4) wherek2Z+7!l;ij(k) is a simple ltering of the input signal u(k) by the Kautz function ;ij, i.e., l;ij(k) =V1X ni=0 ;ij(ni)u(kni); (5) whereV= maxfJ1;:::;Jg. More information on the identi cation approach based on Volterra series can be found in [40, 41]. Details about Kautz functions are given in section 2.2. The reader is also encouraged to see [42, 43, 37]. 2.2. Deterministic Kautz functions The generalized form of the Kautz functions is written as [44] ;2j1(z) =zq (1b2)(1c2) z2+b(c1)zccz2+b(c1)z+ 1 z2+b(c1)zcj1 ; (6) and ;2j(z) = ;2j1(z)zbq 1b2; (7) being the values of bandc, respectively, de ned by b=(Z+Z) 1 +ZZ(8) 6 and c=ZZ; (9) whereZand Zare, respectively, the Kautz poles and its complex conjugate in discrete domain. The discrete poles Zcan be related to the continuous poles Sthrough the equation Z= expS Fs ; (10) whereFsrepresents the sampling frequency. The Kautz poles for each kernel are de ned as S=!j!q 12; (11) where!andare the Kautz parameters and represents the number of the kernel considered. Remembering that, for the linear kernel, the Kautz parame- ters are the natural frequency and damping ratio of the system. In identi cation processes, it is common to use some optimization methodology to nd !and [45]. The use of the Kautz functions as the orthonormal basis is driven by their properties to describe the oscillatory dynamic systems [38]. 2.3. System uncertainties Uncertainties are classi ed in the literature as being of two types, aleatory or epistemic. The rst type is intrinsic to scenarios with variabilities, such as mea- surement noise in the system observations (experimental data) and variabilities of the real system with respect to its nominal con guration (due to geomet- ric imperfections, manufacturing irregularities, environmental conditions, etc.) [4, 5, 6]. These uncertainties can not be eliminated, only better characterized. On the other hand, epistemic uncertainties are due to the lack of knowledge (ignorance) about the system physics. By increasing knowledge about a certain system, these uncertainties can be mitigated [4, 5, 6]. For the sake of simplicity, in this work, only the aleatory uncertainties, also known as data uncertainties, are taken into account. In the context of system identi cation, such uncertainties are materialized in the form of variations in the model parameters. Implicit in this approach is the hypothesis that the Volterra series is capable of producing a reliable representation of the system response. 2.4. The stochastic Volterra series In this work, a stochastic version of the Volterra series is proposed to identify the nonlinear system of interest. The model parameters subjected to uncertain- ties are described as random variables or processes, de ned on the probability space (;;P), where  is a sample space, is a-algebra over , and Pis a probability measure. It is assumed that any random variable 27!Y()2R in this probabilistic setting, with a probability distribution PY(dy) onR, admits a probability density function (PDF) y7!pY(y) with respect to dy. 7 Therefore, the discrete-time Volterra series describes the random nonlinear system output as y(;k) =1X =1N11X n1=0:::N1X n=0H(;n1;:::;n)Y i=1u(kni); (12) whereu(k) is the deterministic input signal, the random process ( ;k)2 Z+7!y(;k) is the aleatory system response and ( ;n1;::;n)2Z7! H(;n1;:::;n) represents the random version of the -order Volterra kernel. The kernels expansion in terms of Kautz functions now is written as H(;n1;:::;n)J1X i1=1:::JX i=1B(;i1;:::;i)Y j=1 ;ij(;nj); (13) whereJ1;:::;Jare the number of Kautz functions used in each orthonormal projection of the Volterra kernels, the random process ( ;i1;:::;i)2Z +7! B(;i1;:::;i) represents the -order random Volterra kernel, expanded in the orthonormal basis, and ( ;nj)2Z+7! ;ij(;nj) represents the random version of the ij-th Kautz lter related with the -order random Volterra kernel. The Kautz functions are modeled as random processes, because their de nition is associated with the dynamic system response y(;k) and depends on the damping ratio and natural frequency, which are subjected to uncertainties. The stochastic version equivalent of the approximation shown in (4) is writ- ten y(;k)1X =1J1X i1=1:::JX i=1B(;i1;:::;i)Y j=1l;ij(;k); (14) where the random process ( ;k)2Z+7!l;ij(;k) is a simple ltering of the deterministic input signal u(k) by the random Kautz function ;ij, i. e., l;ij(;k) =V1X ni=0 ;ij(;ni)u(kni); (15) forV= maxfJ1;:::;Jg. Then, the coecients of the kernels can be estimated, considering Monte Carlo simulations and the least squares method. Figure 1 shows a owchart of the Volterra kernels identi cation considering the random Kautz functions and the Monte Carlo simulations. So, considering each stochastic realization , the matrix can be completed with the input signal ltered l;ij(;k) and the vector ywith the experimental output signal y(;k) = (T)1Ty; (16) where has the terms of the orthonormal kernels B, in each realization . The 8 procedure is repeated until the Monte Carlo convergence is achieved. Figure 1: Description of the Volterra kernels identi cation approach used, based on Monte Carlo simulations. 2.5. The stochastic Kautz parameters In order to create a robust identi cation, it is necessary to provide some type of certi cation, i. e., an envelope of reliability for the nominal values of the model parameters. Such certi cation can be obtained by a stochastic model, where probability distributions, instead of scalar values, are identi ed for the model parameters. In this case, the system response also becomes random, a stochastic process to be precise, allowing the model to predict the variations in the response related to uncertainties. As consequence of the system response variation, the natural frequency and the damping ratio of the equivalent linear system can also vary, becoming random variables 27!n()2R,27!n()2R. As mentioned before, the Kautz parameters for each kernel are related with the systems natural frequency and damping ratio, so they are also considered as random variables 27!()2R,27!()2R. So, the random system response is used to de ne the Kautz poles values. Thus, the Kautz parameters, related with the rst Volterra kernel, can be ap- proximated as natural frequencies and damping ratios of the equivalent linear 9 system. These values can be obtained through modal analysis techniques and using the low amplitude system response, which is approximately linear. As the system linear dynamics vary, it is clear that the modal parameters, and consequently, the Kautz parameters vary too. The de nition of the Kautz pa- rameters related with the high order Volterra kernels is more complicated since they are not exactly natural frequencies and damping ratios. But, it is possible to de ne linear approximations between modal parameters and the high order Kautz parameters 1() =n(); 1() =n(); 2() =p1n(); 2() =p2n(); (17) 3() =p3n(); 3() =p4n(); where the relationships p1;p2;p3andp4are estimated minimizing the error function min J=NX k=1[^y(k)y(k)]2(18) wherey(k) is the deterministic discrete-time system response, ^ y(k) is the de- terministic model response and the vector  contains the coecients of the relationships between modal and Kautz parameters. Once the coecients were de ned, the Kautz parameters variation can be obtained based on the model parameters variation. The optimization procedure does not need to be repeated every time, but only for one representative data. In the simulations performed, in this work, only the response correspondent to the deterministic values of the system parameters was used to optimize the Kautz poles. And, the variation of the Kautz poles was considered through the estimation of the modal parameters in each Monte Carlo realization. 2.6. Damage detection based on Volterra series The stochastic Volterra series can be used as a mathematical model to ap- proximate the system response, rst in the healthy condition to de ne a set of reference models and, then, for each unknown condition (healthy or damaged), to make predictions about the condition of the system. Since the stochastic ver- sion of the Volterra series is able to describe linear and nonlinear systems, the damage detection approach can be applied with linear or nonlinear behavior in reference or damaged conditions. So, the system can be identi ed several times, using Monte Carlo simulations, in the reference condition, to construct a set of Volterra models. In the unknown condition, a new Volterra model can be iden- ti ed and compared with the set of models constructed to detect the possible presence of structure variations related with damages, all this, considering the presence of uncertainties. The stochastic version of the Volterra series proposed allows the damage detection with probability con dence. In this work, two dif- ferent approaches are shown, based on the kernels coecients and on the linear 10 and nonlinear contributions to the total system response, in each structural condition. The approaches are described as follow. 2.6.1. Damage detection based on kernels coecients The rst approach is based on the coecients of the Volterra kernels iden- ti ed. This methodology does not consider the direct in uence of the Kautz functions, because only the coecients expanded in the orthonormal basis are considered. The in uence of the random Kautz functions is taken into ac- count only in the process of model identi cation. Considering Eq. (16) and the Volterra series truncated in the third order term, the reduced order kernels coecients, for each Monte Carlo realization, are allocated as =2 41 2 33 5;1=8 >>>< >>>:B1(1) B1(2) ... B1(J1)9 >>>= >>>;;2=8 >>>>>>>>>>>< >>>>>>>>>>>:B2(1;1) B2(1;2) ... B2(1;J2) B2(J2;1) ... B2(J2;J2)9 >>>>>>>>>>>= >>>>>>>>>>>;;3=8 >>>>>>>>>>>>>>>>< >>>>>>>>>>>>>>>>:B3(1;1;1) B3(1;2;1) ... B3(1;J3;1) B3(1;1;2) ... B3(1;J3;J3) ... B3(J3;J3;J3)9 >>>>>>>>>>>>>>>>= >>>>>>>>>>>>>>>>;;(19) remembering that J1,J2andJ3are the number of Kautz functions used, con- sidering the rst, second and third kernels, respectively. It is expected that the main information about the high order kernels is allocated in the main diago- nal, so it is considered as monitoring parameters only the coecients positioned in this region. Then, the coecients with higher contribution can be used to represent the model identi ed 1=8 >>>< >>>:B1(1) B1(2) ... B1(J1)9 >>>= >>>;;2=8 >>>< >>>:B2(1;1) B2(2;2) ... B2(J2;J2)9 >>>= >>>;;3=8 >>>< >>>:B3(1;1;1) B3(2;2;2) ... B3(J3;J3;J3)9 >>>= >>>;;nl=2 3 :(20) Now, considering the stochastic version of the Volterra series and the Monte Carlo simulations, multiple models are estimated in the reference condition and each vector becomes a random process. The random processes estimated, using Monte Carlo simulations, represents the set of models identi ed in reference condition 1=8 >>>< >>>:B1(;1) B1(;2) ... B1(;J1)9 >>>= >>>;;2=8 >>>< >>>:B2(;1;1) B2(;2;2) ... B2(;J2;J2)9 >>>= >>>;;3=8 >>>< >>>:B3(;1;1;1) B3(;2;2;2) ... B3(;J3;J3;J3)9 >>>= >>>;;nl=2 3 :(21) 11 Additionally, the index mis used to represent the four di erent indicators calculated. So, the general notation mis used from here, with m= 1;2;3; ornl, depending on the number of the kernel considered in the analysis. The damage detection can be summarized as follows 1. Identi cation of the multiple Volterra models in the reference condition, through the stochastic version of the Volterra series; 2. Construction of the indexes based on kernels coecients, in the reference condition ( m); 3. Identi cation of a new Volterra model in an unknown condition; 4. Construction of the indexes based on kernels coecients, in an unknown condition ( m); 5. Comparison between the coecients estimated in unknown and reference conditions, if the new model belongs to the set of reference models, the structure is classi ed as healthy. The comparison between new and reference indexes will be done based on the novelty detection or outliers analysis in multivariate data, using Mahalanobis squared distance, described further. 2.6.2. Damage detection based on kernels contribution With the aim of comparison between the results, a di erent approach is pro- posed, based on the contribution of the Volterra kernels identi ed, to the total response. In this situation, the convolution between the input signal ltered by the Kautz functions and the kernels expanded in the orthonormal basis is considered, giving higher importance to the random Kautz functions. The main advantage of Volterra series approach is the capability of separate the model response in linear and nonlinear contributions, through the kernels estimated y1(k)J1X i1=1B1(i1)li1(k); (22) y2(k)J2X i1=1J2X i2=1B2(i1;i2)li1(k)li2(k); (23) y3(k)J3X i1=1J3X i2=1J3X i3=1B3(i1;i2;i3)li1(k)li2(k)li3(k); (24) ynl(k) =y2(k) +y3(k): (25) wherey1(k),y2(k),y3(k) andynl(k) are the linear, quadratic, cubic and nonlin- ear contributions, respectively. Again, considering the stochastic version of the Volterra series and the Monte Carlo simulations, multiple models are estimated 12 in the reference condition and each vector becomes a random process y1(;k)J1X i1=1B1(;i1)li1(;k); (26) y2(;k)J2X i1=1J2X i2=1B2(;i1;i2)li1(;k)li2(;k); (27) y3(;k)J3X i1=1J3X i2=1J3X i3=1B3(;i1;i2;i3)li1(;k)li2(;k)li3(;k); (28) ynl(;k) =y2(;k) +y3(;k): (29) As mentioned before, the index mis used to represent the ve di erent indexes calculated. So, the general notation ymis used from here, with m= 1;2;3;ornl. The damage detection is very similar to the methodology proposed using kernels coecients and can be summarized as follows 1. Identi cation of the multiple Volterra models in the reference condition, through the stochastic version of the Volterra series; 2. Calculation of the kernels contribution to the total model response, in the reference condition ( ym); 3. Identi cation of the Volterra model in an unknown condition; 4. Calculation of the kernels contribution to the total model response, in an unknown condition ( ym); 5. Comparison between the contributions estimated in unknown and refer- ence condition, if the new model belongs to the set of reference models, the structure is classi ed as healthy. The di erence here is that using the kernels contribution, the Kautz func- tions have an in uence on the process, whereas using only the coecients, they do not. Based on that, the variations in the frequency of oscillation and the damping coecients are captured by the Kautz parameters, this version of the approach is sensitive to these variabilities. The damage detection is also done based on the novelty detection, considering multivariate data. 2.6.3. Novelty detection The indexes proposed based on Volterra models are described as random processes in the reference condition. This methodology combined with Monte Carlo simulations allows the generation of a set of reference models that can be used as monitoring parameters, as described before. So, the damage detection can be done based on the novelty detection or outliers analysis, considering multivariate data. Let us consider the set of indexes in the reference condition (morym), the classical Mahalanobis squared distance can be calculated [15] Dm() = [m(;i)m]T-1[m(;i)m] (30) 13 or Dm() = [ym(;k)ym]T-1[ym(;k)ym] (31) whereDm() is the Mahalanobis squared distance in the reference condition,  is the covariance matrix and (:)is the mean operator. Note that, in the ref- erence condition, Dm() is a random variable, with a correspondent probability density function (PDF). So, its PDF can be empirically estimated considering the Kernel Density Estimation approach in order to establish a threshold value to the distribution [46, 13] ^pDm(dm) =1 NhNX i=1Kdmdm;i h ; (32) where ^pDm(dm) is an approximation of pDm(dm) that is the true density of Dm(),Nis the total number of realizations of the random parameter, dm;i is theith realization of the random variable Dm(),Krepresents the kernel of the transformation, in this work the Gaussian kernel is used, and his the smoothing parameter that controls the width of the Gaussian kernel chosen. The main diculty in the Kernel Density Estimation approach is the choice of the optimal value to the smoothing parameter and, in this work, the cross-validation was used to set the better value to h[47]. The Kernel Density Estimation is used, because the distribution of Dm() is assumed as unknown a priori and the large number of Monte Carlo simulations, provides sucient samples to estimate its density using this approach. Then, with the density estimated ^ pDm(dm) it is possible to establish a thresh- old value for the distribution, which can be used in the damage detection pro- cedure. The threshold value can be de ned as [11] m=fdmsuch thatZ+1 dmpDm(dm) ddm= g; (33) where mrepresents the threshold value and is the sensitivity of the approach or probability of false alarms considered. The better value of to be chosen it depends on the level of security/con dence that each operational application demands. Now, a new model is identi ed in an unknown condition (healthy or damaged). The indexes are computed considering the new model ( unk moryunk m) and compared with the stochastic model (set of reference models), through the Mahalanobis squared distance Dunk m= [unk mm]T-1[unk mm] (34) or Dunk m= [yunk mym]T-1[yunk mym] (35) 14 whereDunk mis the Mahalanobis squared distance in the unknown condition. Finally, the hypothesis test can be applied to determine if the system is in healthy or damaged condition H0:Dunk mm; H1:Dunk m>m;(36) where the null hypothesis H0represents the healthy condition and H1the dam- aged. The methodology used to detect damage is summarized in the owchart showed in g. 2, considering the two proposed approaches. Figure 2: Description of the damage detection approach based on stochastic Volterra series. 3. Modeling a breathing crack in an uncertain nonlinear beam This section describes the studied nonlinear system that can be approxi- mated by a nonlinear single-degree-of-freedom (SDOF) model, a Dung oscil- lator [48]. The deterministic reference system is presented, then, the simulated damage used, a breathing cracked model, is shown. Finally, the upgrade to the stochastic version is done with the necessary information about the uncertainties considered. 3.1. Deterministic system in reference condition In this work, a single-degree-of-freedom (SDOF) model is used to simulate a real system, composed by a clamped-free aluminum beam with a steel mass po- sitioned in its free extremity, to generate a nonlinear interaction with a magnet ( g. 3). The objective of the setup used is to emulate a hardening nonlinear 15 system, to study the performance of the damage detection approach proposed. In all simulations, the input and output signals considered are the force applied by the shaker and the velocity measured near to the free extremity of the beam, respectively. The real system presents nonlinear behavior only for large dis- placements and its dynamic behavior can be well approximated, for a limited frequency range, by its rst mode shape that can be represented as a Dung oscillator [48] mx(t) +c_x(t) +k1x(t) +k2x(t)2+k3x(t)3=U(t); (37) wheremis the system equivalent mass in [kg], cis the damping coecient in [Ns/m],k1is the linear sti ness in [N/m], k2is the quadratic sti ness in [N/m2], k3is the cubic sti ness in [N/m3], andU(t) is the external force in [N]. The displacement, velocity and acceleration in the free extremity of the beam are, respectively, represented by x(t), _x(t) and x(t). The quadratic sti ness was added to emulate the interaction between shaker and beam, to ensure a more realistic description of the system behavior, whereas the cubic sti ness describes the interaction between beam and magnet. Figure 3: Nonlinear system simulated in the reference condition. The deterministic values of the Dung oscillator parameters were estimated and are shown in Table. 1. The estimation of the parameters was performed considering modal analysis and the restoring force surface (RFS) method, con- sidering experimental data measured. The description of the procedure used to approximate the real system considered by the Dung oscillator model it is not shown here because it is not the focus of this work. In order to illustrate the nonlinear behavior of the system in the reference condition (Healthy), some tests were carried out. In all the tests, it was used a sampling frequency of 512 Hz, 2048 samples and it was considered two levels of the input signal (0.1 N - Low level and 1 N - High level) in order to study the linear and nonlinear behavior of the system. The system was excited by a chirp signal varying the excitation frequency from 15 to 30 Hz in 4 seconds. Figure 4 shows the time-frequency diagram of the system response signals, considering the low and high level of the input amplitude. The nonlinear e ect can be observed with the appearance of multiple harmonics (2 and 3 times the fundamental harmonic) when a high level of input is applied. These results 16 Table 1: Dung oscillator parameters. Parameter Deterministic value m[kg] 0.26 c[Ns/m] 1.36 k1[N/m] 5:49103 k2[N/m2] 3:24104 k3[N/m3] 4:68107 con rm that the system has linear behavior to the low level of input signal and nonlinear to the high level of input, in the reference condition. The presence of multiple harmonics is related to the nonlinear sti nesses of the Dung oscillator motion equation. These results show, clearly, the nonlinear behavior of the system for large displacements (condition reached considering 1 N as excitation amplitude). Frequency [Hz] Amplitude [(ms-1)2/Hz] Time [s] 150 130 110 90 70 50 30 10 0.5 1 1.5 2 2.5 3 3.5-18-16-14-12-10-8-6-4 (a) Low input (0.1 N). Frequency [Hz] Amplitude [(ms-1)2/Hz] Time [s] 150 130 110 90 70 50 30 10 0.5 1 1.5 2 2.5 3 3.5-14-12-10-8-6-4-2 (b) High input (1 N). Figure 4: Time-frequency diagram of the response for di erent levels of input amplitude, considering the reference condition (Healthy). 3.2. Damage simulated A breathing crack model The methodology showed in section 2 proposes to detect damages, or struc- tural variations, in nonlinear systems, considering data variation. In this sense, a breathing crack model is used to emulate the presence of a damage in the reference studied system. Figure 5 shows the presence of a crack in the system of interesting. The objective is to detect the presence of the crack in the beam, consider- ing that the system has nonlinear behavior before (Dung oscillator) and after (breathing crack) the damage occurrence, taking into account the data varia- tion. In other words, the nonlinear characteristic of the damage has di erent nature of the nonlinear behavior of the system in the reference condition. A 17 Figure 5: Nonlinear system simulated in the damaged condition. lot of papers, over the years, have shown the approximation of the breathing crack phenomenon by a single-degree-of-freedom (SDOF) model, as a bilinear oscillator [49, 50, 34, 51]. Combining the reference model (Dung oscillator) with the cracked beam model, the system behavior can be approximated by mx(t) +c_x(t) +F(x(t);t) +k2x(t)2+k3x(t)3=U(t); (38) whereF(x(t);t) is the bilinear force F(x(t);t) =( k1x(t); ifx(t)0 k1x(t);ifx(t)<0 where 0< < 1 represents the crack severity. If = 1 the system is in the ref- erence condition (without crack). The breathing crack phenomenon generates a nonlinear behavior with di erent nature of the cubic e ect of the reference model but similar to the quadratic e ect of k2. This situation makes the application of damage detection procedures dicult. The applied approach has to be able to detect the nonlinear damage without confusions with the nonlinear behavior of the reference condition. This is a great improvement compared with other authors like Peng et al. (2007) [52], Surace et al. (2011) [53] and Prawin et al. (2018) [51], for example, that have used the Volterra series to detect damage, considering the linear behavior in reference condition and the nonlinear e ects as a consequence of the crack occurrence. In the presence of the breathing crack, the nonlinear behavior does not de- pend on the amplitude of the response and can be observed for all excitation amplitudes. Figure 6 shows the time-frequency diagram of the system response signals to the low and high level of amplitude excitation, considering the pres- ence of a crack with = 0:9. It can be observed, that to the low and high level of input amplitude, the appearance of quadratic harmonics, mainly of the second order. This behavior brings closer the breathing crack phenomenon by a SDOF model and it is a consequence of the bilinear characteristic of the response. 18 Frequency [Hz] Amplitude [(ms-1)2/Hz] Time [s] 150 130 110 90 70 50 30 10 0.5 1 1.5 2 2.5 3 3.5-16-14-12-10-8-6-4 (a) Low input (0.1 N). Frequency [Hz] Amplitude [(ms-1)2/Hz] Time [s] 150 130 110 90 70 50 30 10 0.5 1 1.5 2 2.5 3 3.5-14-12-10-8-6-4-2 (b) High input (1 N). Figure 6: Time-frequency diagram of the response for di erent levels of input amplitude, considering the crack severity = 0:9. 3.3. Stochastic version of the nonlinear system It is known that mechanical systems and structures are subjected to un- certainties that have to be considered in the processes of models identi cation and damage detection procedures, as it was exempli ed in section 1. So, tak- ing the uncertainties into account is essential to make robust predictions [54]. The data uncertainties considered are related with environmental or operating conditions, temperature and humidity changes, sensor bonding conditions, that can generate variations in the fundamental frequencies and damping ratios of the structures and systems, making dicult the application of damage detec- tion procedures [55]. The main idea is to simulate these variabilities in order to reproduce, in a more realistic way, the diculties involved in the damage detec- tion procedures. In this context, it is not considered the simple noise addition in the measurements, but the variations in the system dynamic behavior. In order to simulate the presence of uncertainties in the system and to test the performance of the stochastic identi cation and damage detection approach proposed in section 2, the linear sti ness ( k1) and damping coecient ( c) of Eq. (38) are considered as random variables 27!k1()2R,27!c()2R. The variation of k1andcre ects directly in the system dynamics, since these parameters in uence the fundamental frequency of oscillation and the damping ratio of the system. So, the Eq. (38) can be rewritten as mx(;t) +c() _x(;t) +F(x(;t);t) +k2x(;t)2+k3x(;t)3=U(t);(39) where the bilinear random force F(x(;t);t) can be described as F(x(;t);t) =( k1()x(;t); ifx(;t)0 k1()x(;t);ifx(;t)<0 and the random processes ( ;t)2R7!x(;t), (;t)2R7!_x(;t), and 19 (;t)2R7!x(;t), represent, respectively, the displacement, velocity and acceleration in the beam free extremity. To construct the probabilistic model, the PDFs of the random parameters were determined through the maximum entropy principle [56, 57]. Considering that the linear sti ness and the damping coecient can not be negative, we assumed the interval (0 ;1) as the support of these random variables. Also, it was considered that the expected value of k1andcare known real numbers k1 andc. For technical reasons, see Soize (2017) [6] for details, we also supposed that the expected value of ln k1and lncare nite. Using these conditions as known information, such as done in Cunha and Sampaio (2015) [58], the principle of maximum entropy says that the probability density function (PDF) of these random variables is given by pZ(z) =1(0;1)1 Z1 2 Z 1 2 Z 1 (1=2 Z)z Z 1 2 Z1 exp z 2 ZZ (40) whereZrepresents the random parameter ( k1orc),Zis the mean value, Zis the dispersion,  indicates the gamma function, and 1(0;1)denotes the indicator function of the interval (0 ;1). This PDF corresponds to a gamma distribution. The stochastic formulation was made assuming, for the lack of better knowledge, that k1andcare independent random variables. To compute the propagation of the uncertainties of the random parameters k1andcthrough the model, the Monte Carlo (MC) method [59] was employed, using 64 samples of each parameter, combined to generate a total of 2048 MC realizations of the deterministic system. The simulations were performed considering k1= 5:49 103[N/m],k1= 0:01,c= 1:36 [Ns/m], c= 0:01. These values of dispersion were chosen to generate satisfactory variation in the system response. Figure 7 shows the probability densities of the random parameters used in the MC simulations. The variation of the random parameters makes the system response varies too and the data variation will be considered in the model identi cation and damage detection procedures. To exemplify the diculty to detect damage in this situation, g. 8 shows the 99% con dence bands of the system Frequency Response Function (FRF), considering 3 di erent structural conditions. It is not possible to distinguish the damaged conditions of the reference state, considering the variation of the response, since it can be observed an overlap of the envelopes. So, the use of modal parameters or the direct system response to detect damage in this situa- tion is not recommended, because only conditions of high severity damage can be detected with probabilistic con dence. The propose is to use the nonlinear behavior of the damage as an indicator of the structural variation. However, the nonlinear behavior in the reference condition has to be considered to avoid the confusion with damage. Next section shows the application of the methodology proposed to detect damages, based on stochastic Volterra series, in the cracked beam model de- scribed, considering the initially nonlinear behavior and data variation. 20 00.10.20.30.40.5 -3 -2 -1 0 1 2 35.40 5.43 5.46 5.49 5.52 5.55 5.58Probability Density Function Normalized ValueNominal Value [kN/m] (a) Linear Sti ness.00.10.20.30.40.5 -3 -2 -1 0 1 2 31.33 1.34 1.35 1.36 1.37 1.38 1.39Probability Density Function Normalized ValueNominal Value [Ns/m] (b) Damping coecient. Figure 7: Distribution of the random parameters k1andc. 10-1100 16 18 20 22 24 26 28 30Amplitude [(ms-1)/Hz] Frequency [Hz] (a) Low level of input (0.1 N).10-1100 16 18 20 22 24 26 28 30Amplitude [(ms-1)/Hz] Frequency [Hz] (b) High level of input (1 N). Figure 8: FRFs considering di erent structural conditions and 2 levels of input amplitude. B represents = 1 (reference condition), B = 0:98 (damaged condition) and B = 0:96 (damaged condition). 4. Damage detection procedure This section shows the principal results obtained with the application of the damage detection approach proposed, based on the stochastic Volterra se- ries, in the nonlinear system described in section 3, considering the nonlinear behavior of the system in the reference condition and the data variation. In all simulations, it was used a sampling frequency of 512 Hz and 2048 sam- ples. The Monte Carlo simulations were performed considering 2048 realiza- tions, to ensure the method convergence. The response of the system described in section 3 was simulated considering di erent structural conditions, i. e., = 1:00;0:98;0:96;0:94;0:92;0:90;0:88;0:86 (2 reference and 7 damage con- ditions) and the simulated data variation. The simulations regarding the refer- 21 ence condition were repeated twice, 2048 realizations were used as training data and 2048 as test data. To better emulate a real application, it was added to the system response, in all situations simulated, a Gaussian noise in order to have a signal to noise ratio (SNR) of 30 dB. 4.1. Stochastic Volterra models identi cation The identi cation of the stochastic Volterra model was done considering 2048 realizations, for each structural condition, obtained through the simulations per- formed. Only the rst three kernels were used to identify the system, once its nonlinear characteristic in the healthy and damaged condition can be approx- imated by the quadratic and cubic term. First of all, the Kautz parameters were parametrized in terms of the natural frequency ( n()) and the damping ratio (()) of the equivalent linear system. As the uncertainties considered in this work in uence on the modal parameters, they vary in each realization. So, they were estimated to each response realization and the optimization process described in section 2.5 was used to de ne the Kautz parameters, taking into account the relations described in Eqs. (17). Table 2 shows the values of the factorsp1,p2,p3andp4found for the high order kernels, considering the dif- ferent damage severities. Changes are observed in the position of the poles as the crack increases, mainly related with the second order kernel, because of the quadratic e ect of the crack. The number of functions were chosen as done in Shiki et al. (2017) [37], and de ned as J1= 2,J2= 4 andJ3= 6. Table 2: Relation factors between the Kautz and modal parameters, considering the high order Volterra kernels. Crack severity Second kernel Third kernel ( ) p1p2p3p4 1.00 1.11 2.7 1.06 1.1 0.98 1.07 2.4 1.06 1.1 0.96 1.07 2.3 1.06 1.1 0.94 1.06 2.2 1.06 1.1 0.92 1.06 2.1 1.06 1.1 0.90 1.05 2.0 1.06 1.0 0.86 1.04 1.8 1.05 1.0 With the Kautz poles de ned, the random Kautz functions can be obtained and used in the models' identi cation process. The identi cation of the Volterra models was performed considering a chirp signal as input applied to the system, varying the excitation frequency from 15 to 30 Hz ( rst mode region). The kernels identi cation was carried out in two steps, rst, the linear kernel was estimated considering low level of input (0.1 N), and then, the high order kernels were estimated considering a high level of input (1 N), as done in Shiki et al. (2017) [37]. The models were estimated via MC simulation, in a total of 2048 realizations for each structural condition. 22 The convergence of the MC method is measured with the aid of a function that depends on the kernel vector or its main diagonal (high order kernels) represented in the time domain, de ned by conv(N) =vuut1 NNX n=1Ztf t=t0jjh(n;t)jj2dt; (41) whereNis the number of MC simulations, jj:jjdenotes the standard Euclidean norm and h(n;t) represents the n-th realization of the random rst kernel or main diagonal of high order kernels, in the time domain. For further details, the reader is encouraged to see Soize (2005) [60]. The criterion was applied to the rst, second and third kernels. Figure 9 shows the results obtained, considering the Volterra kernels identi ed in the reference condition. It can be seen that the MC convergence is achieved with the number of samples used. 2.2752.2802.2852.2902.2952.3002.305 0 500 1000 1500 2000conv (N)103 N (a) First kernel.1.01.52.02.53.03.54.04.55.05.56.0 0 500 1000 1500 2000conv (N)107 N (b) Diagonal of second kernel. 00.51.01.52.02.53.03.5 0 500 1000 1500 2000conv (N)106 N (c) Main diagonal of third kernel. Figure 9: MC convergence test applied to the Volterra kernels identi ed in the reference condition ( = 1:00). 23 4.2. Stochastic Volterra models validation Before the application of the Volterra series models in the damage detection process, it is important to verify if the Volterra kernels describe the dynamical behavior of the system with certain statistical con dence. Thus, the same input used to estimate the Volterra kernels was applied to the model. Figure 10 shows the comparison between the 99% con dence bands of the stochastic Volterra model and new simulated data, considering high level of input (1 N) and two structural conditions, reference ( = 1:00) and severe damage ( = 0:86). It can be seen that the model is able to describe the system behavior and to predict the response, even in the presence of nonlinear behavior and data variation, in both structural conditions. -0.6-0.4-0.200.20.40.6 2 2.2 2.4 2.6 2.8 3 3.2 3.4Velocity [m/s] Time [s] (a) Complete Signal ( = 1:00).-0.8-0.6-0.4-0.200.20.40.60.8 2.80 2.82 2.84 2.86 2.88 2.90Velocity [m/s] Time [s] (b) Zoom ( = 1:00). -0.6-0.4-0.200.20.40.6 2 2.2 2.4 2.6 2.8 3 3.2 3.4Velocity [m/s] Time [s] (c) Complete Signal ( = 0:86).-0.8-0.6-0.4-0.200.20.40.60.8 2.50 2.52 2.54 2.56 2.58 2.60Velocity [m/s] Time [s] (d) Zoom ( = 0:86). Figure 10: Stochastic Volterra model response in comparison with new simulated data, con- sidering a high level amplitude chirp input (1 N) and two structural conditions. The gray box represents the 99% con dence bands, { { represents the mean and { { the simulated data. Then, a signal with di erent nature was applied to the system to validate the stochastic Volterra model in the frequency domain. A single tone sine was applied with a high level of input (1 N) and frequency close to the equivalent linear system natural frequency ( 23 Hz), in order to amplify the nonlinear behavior. Figure 11 shows the results obtained in the frequency domain, to help 24 the visualization of the harmonic components, with 99% of statistical con dence in comparison with new simulated data. Again, two structural conditions are shown, healthy ( = 1:00) and severe damage ( = 0:86). It can be seen that the model predicts the system behavior in healthy ( g. 11a) and damaged ( g. 11b) conditions. All frequency components present in the signal can be described by the stochastic Volterra model. The large dispersion observed is caused by the number of Kautz functions used on the description of the third order Volterra kernel, that ampli es the propagation of the uncertainties into the model, but a lower number of functions does not do the model able to describe the nonlinear behavior of the studied system. With the stochastic Volterra model validated, the stochastic response can be used to detect damage in the system as shown in the next section. 10-810-710-610-510-410-310-210-1100 10 20 30 40 50 60 70 80Amplitude [(ms-1)2/Hz] Frequency [Hz] (a) Healthy condition ( = 1:00).10-810-710-610-510-410-310-210-1100 10 20 30 40 50 60 70 80Amplitude [(ms-1)2/Hz] Frequency [Hz] (b) Damaged condition ( = 0:86). Figure 11: Stochastic Volterra model response in comparison with new simulated data, in frequency domain, considering a single tone sine input and two structural conditions. The gray box represents the 99% con dence bands, { { represents the mean and { { the simulated data. The main advantage in the use of the approach based on Volterra series is the capability to separate linear and nonlinear contributions in the total response, through the Volterra kernels. In this work, the idea is to use this capability to lter and compare the nonlinear contributions to the total system response before and after the damage occurrence. In order to exemplify this information, g. 12 shows the contributions of the rst, second and third kernels with 99% of statistical con dence, considering a high level of input signal (1 N) and two structural conditions, healthy ( = 1:00) and severe damage ( = 0:86). All the contributions, linear, quadratic and cubic, change with the occurrence of the severe damage because the natural frequency and, consequently, the Kautz poles are in uenced by the crack behavior. This result shows that the kernels contributions, or kernels coecients, can be used in the process of damage de- tection. The variation of the quadratic contribution is bigger, because of the nonlinear nature of the damage. Finally, in the situation of severe damage, the di erence between the response in healthy and damaged condition can be visu- 25 ally observed, but at the beginning of the crack propagation, the di erentiation is more dicult, mainly in the presence of uncertainties. Therefore, an index has to be used to detect the structural variations. It is expected that the index related with the second kernel will be more sensitive to the presence of the crack. -0.8-0.6-0.4-0.200.20.40.60.8 1.5 2.0 2.5 3.0 3.5Velocity [m/s] Time [s] (a) Linear contribution.-0.020-0.015-0.010-0.00500.0050.0100.0150.020 1.5 2 2.5 3 3.5Velocity [m/s] Time [s] (b) Quadratic contribution. -0.8-0.6-0.4-0.200.20.40.60.8 1.5 2.0 2.5 3.0 3.5Velocity [m/s] Time [s] (c) Cubic contribution. Figure 12: Con dence bands of the kernels contribution to the total system response obtained through the Volterra model with 99% of con dence, considering high level of input (1 N). B represents the reference condition ( = 1:00) and B the damaged condition ( = 0:86). 4.3. Damage detection This section exempli es the two di erent metrics proposed to detect damage in nonlinear systems, using the stochastic version of the Volterra series. As men- tioned before, the tests were performed considering = 1:00;0:98;0:96;0:94, 0:92;0:90;0:88;0:86 (2 reference and 7 damage conditions) and the simulated data variation, considering the uncertainties in linear sti ness and damping co- ecient. The rst 2048 realizations in the reference condition were used as training data and the others as test data. 4.3.1. Use of kernels coecients The approach proposed, based on Volterra kernels coecients, showed in sec- tion 2.6.1 was applied, combined with the novelty detection showed in section 26 2.6.3. The multiple models identi ed in all structural conditions were considered in the calculation of the indexes. Four di erent indexes were calculated ( 1,2, 3andnl), to compare the results. Figure 13 shows the evolution of the Ma- halanobis squared distance applied to the indexes, with the crack propagation. The boxplot method is used to clarify the visualization of the results. The rst two boxes represent the reference condition, the data used to train the model and the data used to test the model. It is observed that the linear index is not able to detect the crack evolution, considering only the kernels coecients as an indicator. The cubic index has a small increase with the crack evolution, whereas the quadratic index has a larger increase. This occurs because of the quadratic nature of the crack behavior. How expected, the better performance is obtained through the use of the nonlinear index, considering both the coef- cients of the second and third kernels. The dispersion of the values is small, because the uncertainties in the modal parameters have lower in uence in the estimation of the higher order Volterra kernels coecients. The dispersion ob- served in the quadratic index is related to the noise added because the quadratic component has the level of amplitude close to the noise. 101100101102 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Mahalanobis squared distance Crack severity ( ) (a) Linear index ( 1).100101102103 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Mahalanobis squared distance Crack severity ( ) (b) Quadratic index ( 2). 100101102103 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Mahalanobis squared distance Crack severity ( ) (c) Cubic index ( 3).100101102103 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Mahalanobis squared distance Crack severity ( ) (d) Nonlinear index ( nl). Figure 13: Boxplot of the Mahalanobis squared distance applied to the four indexes calculated, for di erent levels of crack severity. With the Mahalanobis squared distance of the indexes calculated in the 27 reference condition, and considering their empirical distribution obtained, it is possible to de ne threshold values, taking into account a probability of false alarms ( ) required in each application. In this work, three di erent values were considered, = 0:005, 0:01 and 0:02. The results of the hypothesis test applied are shown in g. 14, for the three values of used. The performance of the linear index it is not satisfactory, with a low percentage of detection in the initial propagation of the crack. The cubic index can detect the crack depending on the severity of the damage, i. e., for close to 1.00 the index fails, but for 0:86, all damaged conditions are detected. Finally, the quadratic and nonlinear (sum of quadratic and cubic coecients) indexes have similar behavior, detecting the damage even on the initial conditions of crack propagation. Remembering that it was considered the data variation, related with the presence of uncertainties, in the process of kernels identi cation, but, nevertheless, the nonlinear index was able to make di erence between the variations related with changes in the natural frequency and damping ratio and the variations related with the crack behavior. This is an important result, showing that the approach can detect damage in uncertain nonlinear systems with probabilistic con dence. 020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (a) Linear index ( 1).020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (b) Quadratic index ( 2). 020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (c) Cubic index ( 3).020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (d) Nonlinear index ( nl). Figure 14: Percentage of damage detection obtained using the kernels coecients, considering di erent crack severities and probability of false alarms used. B represents = 0:02,B represents = 0:01 and B represents = 0:005. 28 With the aim of study the performance of the approach proposed, the re- ceiver operating characteristics (ROC) curve was computed [2]. This curve presents the relation between true detections and false alarms for di erent values. Figure 15 shows the results obtained for the 4 indexes. The linear index has the worst performance and the cubic index an intermediate perfor- mance. The quadratic and nonlinear indexes have similar performance, with a higher rate of detection without expressive false alarms rate. These results show that the nonlinear index proposed is able to detect the initial propagation of the crack even in the presence of uncertainties, simulated by variations in the natural frequency and damping ratio of the system. This high performance is possible because of the nonlinear nature of the damage simulated that has a large in uence in the coecients of the second kernel estimated. The use of the Volterra series in this situation improve the results in the damage detection process, representing a real contribution when the data variation related to the uncertainties is considered. 00.20.40.60.81.0 0 0.2 0.4 0.6 0.8 1.0True detection rate False alarm rate Figure 15: Receiver operating characteristics (ROC) curve, obtained using the kernels coe- cients. { { represents the linear index, { x { the quadratic index, { { the cubic index and {{ the nonlinear index. 4.3.2. Use of kernels contribution Now, the approach proposed based on Volterra kernels contribution showed in section 2.6.2 was applied, combined with the novelty detection showed in sec- tion 2.6.3. This methodology uses the main advantage of the use of the Volterra series, that is the capability of separate the kernels contribution to the total re- sponse. The di erence compared with the approach presented in the last section is the in uence of the Kautz functions in the process, that introduces the in u- ence of the natural frequency and damping ratio in the approach, through the Kautz poles dependence. The increase of the crack makes the natural frequency 29 of the system change, and this variation can be detected using this approach, but, on the other hand, it is expected that the uncertainties have more in uence using this approach, doing the indexes vary more. Again, the multiple models identi ed in all system conditions (healthy and damaged) were considered in the calculation of the contributions. The input used here is the same chirp signal used in the estimation of the kernels, con- sidering a high level of amplitude (1 N), to maximizes the nonlinear in uence in the response. Four di erent indexes were calculated ( y1,y2,y3andynl), to compare the results. Figure 16 shows the evolution of the Mahalanobis squared distance applied to the kernels contribution, with the crack progression. It is clear that using this approach the linear and cubic indexes have better perfor- mance to detect the presence of the crack compared with the use of kernels coecients. The increase of the indexes is related to the variation of the natu- ral frequency with the progression of the damage. Again, the better results are obtained through the use of the quadratic or nonlinear indexes, because of the nature of the damage simulated. The use of the kernels contribution increases the performance of the linear and cubic indexes, but the performance of the quadratic and nonlinear indexes continued to be better. With the Mahalanobis squared distance of the kernels contribution calcu- lated in the reference condition, and considering their empirical distribution obtained, it is possible to de ne threshold values, taking into account the prob- ability of false alarms ( ) required in each application. Again, three di erent values were considered, = 0:005, 0:01 and 0:02. The results of the hypothesis test applied are shown in g. 17, for the three values of used. The perfor- mance of the linear and cubic indexes was improved with the use of kernels contribution, with a higher percentage of detection of the crack for 0:92. The quadratic and nonlinear (sum of quadratic and cubic contributions) indexes were able to detect the damage even on initial conditions of crack propagation. The variation of the frequency and damping is considered in the Kautz poles estimation, improving the capability of the approach to detect the frequency variation related with the crack behavior, but confusing this variation with the changes that are a consequence of the presence of uncertainties. Finally, the receiver operating characteristics (ROC) curve was computed again [2]. Figure 18 shows the results obtained in all situations, considering linear and nonlinear kernels. The linear and cubic contributions have the worst performance again because these components of the response have a slight vari- ation for values of crack close to 1.00. The quadratic and nonlinear indexes have close performance, with a higher rate of detection without expressive false alarms rate. These results show that the analysis of the nonlinear contribution to the total response is able to detect the initial propagation of the crack even in the presence of the uncertainties. This higher performance is possible because of the nonlinear nature of the damage simulated, that has a large in uence on the nonlinear dynamics of the system. Again, the use of the Volterra series in this situation represents an advantage in the damage detection process, with the improvement in the results in comparison with the linear approach represented by the rst kernel contribution. 30 100102104106108101010121014 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96Mahalanobis squared distance Crack severity ( ) (a) Linear contribution ( y1).100102104106108101010121014 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96Mahalanobis squared distance Crack severity ( ) (b) Quadratic contribution ( y2). 100102104106108101010121014 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96Mahalanobis squared distance Crack severity ( ) (c) Cubic contribution ( y3).100102104106108101010121014 1.00 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96Mahalanobis squared distance Crack severity ( ) (d) Nonlinear contribution ( ynl=y2+ y3). Figure 16: Boxplot of Mahalanobis squared distance of the kernels contribution calculated for di erent levels of crack severity. 4.3.3. Comparison between the use of kernels coecients and kernels contribu- tion Comparing the results that have been shown in gs. 13 and 16, it can be noticed that the indexes variation is bigger using the kernels contribution, be- cause the data variation induces the Kautz poles to vary, making the use of kernels contribution more sensitive to the presence of uncertainties than the use of the kernels coecients. On the other hand, the capability of the Kautz func- tions to describe the variation of the frequency of oscillation related with the crack presence makes the approach more sensitive to the presence of the dam- age, mainly considering the linear component of the response. So, damages that induce linear variations in the response, related with frequencies of oscillation or damping ratios, can be better detected using the kernels contribution. Addi- tionally, comparing the gs. 15 and 18 one can observe that the performance of the linear and cubic indexes improves with the use of the kernels contribution, with a higher rate of true detection, but in some sceneries where the presence of the uncertainties can have higher in uence in the system response, the confusion between the data variation and the presence of damage must be greater. Now, when the damage induces the system to present nonlinear behavior, 31 020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (a) Linear contribution ( y1).020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (b) Quadratic contribution ( y2). 020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (c) Cubic contribution ( y3).020406080100 1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86Percentage of detection [%] Crack severity ( ) (d) Nonlinear contribution ( ynl=y2+ y3). Figure 17: Percentage of damage detection obtained using the kernels contribution, consider- ing di erent crack severities and probability of false alarms used. B represents = 0:02,B represents = 0:01 and B represents = 0:005. like the breathing crack phenomenon studied, the coecients of the higher order kernels showed to be sensitive to the presence of the damage. In this situation, the approach based on the kernels coecients is more e ective to make di er- ence between the uncertainties and the damages. Therefore, both approaches have shown interesting performance and di erent characteristics, related to the sensitivity to the presence of the damage and the uncertainties. So, the decision about the use of one or other approach it will depend on the problem consid- ered, given that, when the damage induces nonlinear behavior to the system response, the level of noise it is not so high and the uncertainties are related with linear components, like the example simulated in this paper, the use of the kernels coecients it is interesting. 5. Final remarks The paper has presented an approach to detect damages in initially non- linear systems, considering data variation caused by the uncertainties. The approach is based on a stochastic version of the Volterra series expanded using 32 00.20.40.60.81.0 0 0.2 0.4 0.6 0.8 1.0True detection rate False alarm rate Figure 18: Receiver operating characteristics (ROC) curve, obtained using the kernels contri- bution. { { represents the linear index, { x { the quadratic index, { { the cubic index and { { the nonlinear index. random Kautz functions. Two di erent methodologies were presented, con- sidering the Volterra kernels coecients and the Volterra kernels contribution. The nal decision about the state of the system (healthy or damaged) was per- formed considering the Mahalanobis squared distance and a threshold value, established from the model estimated in the reference condition. Linear and nonlinear analysis were compared, using the capability of the Volterra series to separate the model's contribution. About the results, the two approaches pro- posed have presented similar performance. As expected, the linear analysis fails when the damage has low severity, because of the in uence of the uncertainties in the system response. The nonlinear analysis uses the nonlinear dynamics of the damage to detect it. It is important to highlight that the nonlinear behav- ior is considered, in this work, even in the reference condition, which combined with the presence of uncertainties, makes dicult the damage detection process. But, the nonlinear indexes have shown to be able to di er the initial nonlinear behavior and data uncertainties of the crack behavior, with satisfactory perfor- mance. As future work, the authors intend to apply the procedure considering an experimental setup composed by a nonlinear beam subject to the presence of a breathing crack. Acknowledgments The authors are thankful for the nancial support provided from of S~ ao Paulo Research Foundation (FAPESP), grant number 2012/09135-3 and 2015/25676- 2, from Research Support Foundation of the State of Rio de Janeiro (FAPERJ), grant number E-26/010.002178/2015, and CNPq grant number 307520/2016-1. 33