Luis Gustavo

Giacon Villani

Samuel Da Silva
email: samuel.silva13@unesp.br

Americo Cunha
email: americo@ime.uerj.br

Luis G G Villani
email: luis.villani@unesp.br

Damage detection in uncertain nonlinear systems based on stochastic Volterra series

Keywords: Damage detection, uncertainties quantification, stochastic Volterra series, nonlinear dynamics

Introduction

The process of damage detection in mechanical, aerospace and civil systems and structures, based on damage features and statistical analysis of them, is commonly called as Structural Health Monitoring (SHM) [START_REF] Farrar | An introduction to structural health monitoring[END_REF]. The study and application of SHM techniques are motivated by the potential life-safety and economic impact of their implementation, which contributes to the development of new approaches [START_REF] Farrar | Structural health monitoring: a machine learning perspective[END_REF]. As it is widely known, the damage-sensitive features can be sensitive to confounding effects such as environmental or operating conditions, temperature and humidity changes, sensor bonding conditions, and others [START_REF] Worden | On switching response surface models, with applications to the structural health monitoring of bridges[END_REF]. 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 [START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamics -Recent theoretical advances[END_REF][START_REF] Soize | Stochastic models of uncertainties in computational mechanics[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. These uncertainties can complicate the damage detection problem and have to be considered in the SHM applications, reducing the occurrence of false alarms [START_REF] Mao | Statistical modeling of frequency response function estimation for uncertainty quantification[END_REF].

To overcome this problem, the damage detection can be done based on statistical procedures, taking into account confidence 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) classification, in this case all states are known and the condition of the structure or system is classified into groups of different damage stages or healthy condition; (iii) regression algorithms, in which some variable or process is monitored with confidence limits. The choice of the more interesting approach to be used it depends on the specific problem and on the level of knowledge about the damage (detection, location, assessment, prediction) desired [START_REF] Worden | Structural Health Monitoring: from Structures to Systems-of-Systems[END_REF]. In this sense, a lot of methods can be used, such as linear and nonlinear Principal Component Analysis (PCA -NPCA) [START_REF] Pavlopoulou | Novelty detection and dimension reduction via guided ultrasonic waves: Damage monitoring of scarf repairs in composite laminates[END_REF] , Extreme Value Statistics (EVS) [START_REF] Sohn | Structural damage classification using extreme value statistics[END_REF], Peaks Over Threshold (POT) [START_REF] Rébillat | Peaks Over Threshold-based detector design for structural health monitoring: Application to aerospace structures[END_REF], machine learning algorithms [START_REF] Santos | Machine learning algorithms for damage detection: Kernel-based approaches[END_REF], neural network [START_REF] Worden | Experimental validation of a structural health monitoring methodology: Part i. novelty detection on a laboratory structure[END_REF], Bayesian approaches [START_REF] Figueiredo | A bayesian approach based on a Markov-Chain Monte Carlo method for damage detection under unknown sources of variability[END_REF], Mahalanobis distance [START_REF] Worden | Damage detection using outlier analysis[END_REF][START_REF] Figueiredo | Autoregressive modeling with state-space embedding vectors for damage detection under operational variability[END_REF], and others. In all these applications, the final objective of the methods used are the same, differentiate the uncertainties and variabilities of the damages.

On the other hand, many existing SHM features use the nonlinear characteristics of damages to detect them. In this context, cracks [START_REF] Lim | Reference-free fatigue crack detection using nonlinear ultrasonic modulation under various temperature and loading conditions[END_REF], delamination [START_REF] Mandal | Health monitoring of stiffened metallic plates using nonlinear wave interaction and embedded pzt transducers[END_REF] or rubbing in rotor systems [START_REF] Zeng | Normalized complex Teager energy operator demodulation method and its application to fault diagnosis in a rubbing rotor system[END_REF], induce linear systems to exhibit nonlinear phenomena 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 [START_REF] Worden | A review of nonlinear dynamics applications to structural health monitoring[END_REF]. Thus, the crack detection in beams can be performed based on the nonlinear behavior induced by this type of damage. Andreaeus et al. [START_REF] Andreaus | Experimental damage detection of cracked beams by using nonlinear characteristics of forced response[END_REF] 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ébillat et al. [START_REF] Rébillat | Nonlinear structural damage detection based on cascade of Hammerstein models[END_REF] have used cascade of Hammerstein 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. [START_REF] Lim | Reference-free fatigue crack detection using nonlinear ultrasonic modulation under various temperature and loading conditions[END_REF] have shown a reference-free fatigue crack detection technique, based on nonlinear ultrasonic modulation. In this last work, a sequential outlier analysis was used to identify the crack presence without referring 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 effects 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 [START_REF] Bornn | Damage detection in initially nonlinear systems[END_REF]. 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 [START_REF] Virgin | Introduction to Experimental Nonlinear Dynamics[END_REF][START_REF] Noël | Nonlinear system identification in structural dynamics: 10 more years of progress[END_REF]. 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 [START_REF] Kerschen | Nonlinear system identification in structural dynamics: current status and future directions[END_REF][START_REF] Worden | Nonlinear system identification in structural dynamics: A short (and biased) history[END_REF], restoring force surface methods (RFS) [START_REF] Masri | A nonparametric identification technique for nonlinear dynamic problems[END_REF][START_REF] Noel | Application of the restoring force surface method to a real-life spacecraft structure[END_REF], harmonic balance method [START_REF] Tang | Using nonlinear jumps to estimate cubic stiffness nonlinearity: An experimental study[END_REF] and others. Unfortunately, these nonlinearities have many individualities, in a way that it is difficult to obtain a general model that describes all the structures of interest. Among the different methods for nonlinearities identification, the Volterra series stands out, because it is a generalization of the linear convolution concept, allowing the separation of the system response in linear and nonlinear components [START_REF] Schetzen | The Volterra and Wiener Theories of Nonlinear Systems[END_REF]. The main procedure to estimate the Volterra kernels is the Harmonic Probing method, that was extensively used in system identification problems [START_REF] Chatterjee | Parameter estimation of Duffing oscillator using Volterra series and multi-tone excitation[END_REF][START_REF] Scussel | The harmonic probing method for output-only nonlinear mechanical systems[END_REF] and damage characterization [START_REF] Chatterjee | Structural damage assessment in a cantilever beam with a breathing crack using higher order frequency response functions[END_REF]. 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. [START_REF] Tang | Fault diagnosis approach based on volterra models[END_REF] [START_REF] Tang | Fault diagnosis approach based on volterra models[END_REF] have used the classical version of the Volterra series approach to detect structural variations in a rotor, using the vibration signals of different points as input/output signals. [START_REF] Tang | Fault diagnosis approach based on volterra models[END_REF] have used the classical version of the Volterra series approach to detect structural variations in a rotor, using the vibration signals of different points as input/output signals. The Volterra kernels coefficients and the Principal Component Analysis (PCA) were used as damages features. In this study, the use of the first 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 identified has difficulted the approach implementation using an optimization process. Similarly, [START_REF] Xia | A novel identification method of Volterra series in rotor-bearing system for fault diagnosis[END_REF] [START_REF] Xia | A novel identification method of Volterra series in rotor-bearing system for fault diagnosis[END_REF] 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 different types of structural variation. In both works, the Volterra series coefficients 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 confidence.

Additionally, Shiki et al. (2017) [START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF] have used this Volterra series formulation, based on input/output signals, to detect structural variations related with damage, in a beam with nonlinear behavior, even in the reference condition. 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 environmental 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 first methodology considers only the Volterra kernels coefficients and the second considers the kernels contribution. The approaches are applied in a simulated uncertain nonlinear system that describes the real behavior of a nonlinear beam, and, the damage is simulated by the propagation of a breathing crack, described as a bilinear stiffness oscillator. The results obtained show that the methodologies proposed are able to differentiate 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-age detection through the approach, even in the presence of uncertainties, with probability confidence.

The present paper is organized as follows. Section 2 reviews the fundamental aspects of system identification 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 uncertainties of different 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.

The stochastic Volterra series approach for damage detection

This section describes the proposed methodology for modeling nonlinear systems assuming data uncertainties. The approach is based on a stochastic version of the Volterra series, expanded using Kautz functions as the orthonormal basis. 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. [START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF]. 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.

Deterministic Volterra series

Consider a discrete-time causal nonlinear system with a single output k ∈ Z + → y(k) caused by a single input k ∈ Z + → u(k), with Z + 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) = ∞ η=1 N1-1 n1=0 . . . Nη-1 nη=0 H η (n 1 , . . . , n η ) η i=1 u(k -n i ), (1)
where (n 1 , .., n η ) ∈ Z η + → H η (n 1 , . . . , n η ) represents the η-order Volterra kernel. This formalism is very versatile, in a way that it allows one to represent different types of nonlinear systems using the convolution concept [START_REF] Schetzen | The Volterra and Wiener Theories of Nonlinear Systems[END_REF]. Also, it allows one to split the system output into a sum of linear and nonlinear contributions

y(k) = y 1 (k) linear + y 2 (k) + y 3 (k) + • • • + y η (k) nonlinear .
(

) 2
The main drawback of the approach is the difficulty of convergence of the series when a large number of terms N 1 , ..., N η are used. In order to reduce the number 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 [START_REF] Kautz | Transient synthesis in the time domain, Circuit Theory[END_REF][START_REF] Heuberger | Modelling and identification with rational orthogonal basis functions[END_REF]. In this way, one can write

H η (n 1 , ... , n η ) ≈ J1 i1=1 ... Jη iη=1 B η (i 1 , ... , i η ) η j=1 ψ η,ij (n j ) , (3)
where J 1 , . . . , J η are the number of Kautz functions used in each orthonormal projections of the Volterra kernels, (i 1 , . . . , i η ) ∈ Z η + → B η (i 1 , . . . , i η ) represents the η-order Volterra kernel, expanded in the orthonormal basis, and n j ∈ Z + → ψ η,ij (n j ) represents the i j -th Kautz filter.

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 filtered by the Kautz functions, i.e.,

y(k) ≈ ∞ η=1 J1 i1=1 . . . Jη iη=1 B η (i 1 , . . . , i η ) η j=1 l η,ij (k) , (4)
where k ∈ Z + → l η,ij (k) is a simple filtering of the input signal u(k) by the Kautz function ψ η,ij , i.e.,

l η,ij (k) = V -1 ni=0 ψ η,ij (n i )u(k -n i ) , (5)
where V = max{J 1 , . . . , J η }. More information on the identification approach based on Volterra series can be found in [START_REF] Silva | Non-parametric identification of mechanical systems by Kautz filter with multiple poles[END_REF][START_REF] Silva | Non-linear model updating of a three-dimensional portal frame based on Wiener series[END_REF]. Details about Kautz functions are given in section 2.2. The reader is also encouraged to see [START_REF] Wahlberg | On approximation of stable linear dynamical systems using Laguerre and Kautz functions[END_REF][START_REF] Oliveira | An introduction to models based on Laguerre, Kautz and other related orthonormal functions -part II: non-linear models[END_REF][START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF].

Deterministic Kautz functions

The generalized form of the Kautz functions is written as [START_REF] Wahlberg | System identification using Kautz models[END_REF] Ψ

η,2j-1 (z) = z (1 -b 2 η )(1 -c 2 η ) z 2 + b η (c η -1)z -c η -c η z 2 + b η (c η -1)z + 1 z 2 + b η (c η -1)z -c η j-1 , (6)
and

Ψ η,2j (z) = Ψ η,2j-1 (z) z -b η 1 -b 2 η , (7)
being the values of b η and c η , respectively, defined by

b η = (Z η + Zη ) 1 + Z η Zη (8)
and

c η = -Z η Zη , (9)
where Z and Z are, 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 η F s , (10)
where F s represents the sampling frequency. The Kautz poles for each kernel are defined as

S η = -ξ η ω η ± jω η 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 parameters are the natural frequency and damping ratio of the system. In identification processes, it is common to use some optimization methodology to find ω η and ξ η [START_REF] Da Rosa | Choice of free parameters in expansions of discrete-time Volterra models using Kautz functions[END_REF]. The use of the Kautz functions as the orthonormal basis is driven by their properties to describe the oscillatory dynamic systems [START_REF] Kautz | Transient synthesis in the time domain, Circuit Theory[END_REF].

System uncertainties

Uncertainties are classified in the literature as being of two types, aleatory or epistemic. The first type is intrinsic to scenarios with variabilities, such as measurement noise in the system observations (experimental data) and variabilities of the real system with respect to its nominal configuration (due to geometric imperfections, manufacturing irregularities, environmental conditions, etc.) [START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamics -Recent theoretical advances[END_REF][START_REF] Soize | Stochastic models of uncertainties in computational mechanics[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. 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 [START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamics -Recent theoretical advances[END_REF][START_REF] Soize | Stochastic models of uncertainties in computational mechanics[END_REF][START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF].

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 identification, 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.

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 uncertainties are described as random variables or processes, defined on the probability space (Θ, Σ, P), where Θ is a sample space, Σ is a σ-algebra over Θ, and P is a probability measure. It is assumed that any random variable θ ∈ Θ → Y(θ) ∈ R in this probabilistic setting, with a probability distribution P Y (dy) on R, admits a probability density function (PDF) y → p Y (y) with respect to dy. Therefore, the discrete-time Volterra series describes the random nonlinear system output as

y(θ, k) = ∞ η=1 N1-1 n1=0 . . . Nη-1 nη=0 H η (θ, n 1 , . . . , n η ) η i=1 u(k -n i ) , (12)
where u(k) is the deterministic input signal, the random process (θ, k) ∈ Θ × Z + → y(θ, k) is the aleatory system response and (θ, n 1 , .., n η ) ∈ Θ × Z η → H η (θ, n 1 , . . . , n η ) represents the random version of the η-order Volterra kernel. The kernels expansion in terms of Kautz functions now is written as

H η (θ, n 1 , ... , n η ) ≈ J1 i1=1 ... Jη iη=1 B η (θ, i 1 , ... , i η ) η j=1 ψ η,ij (θ, n j ) , (13)
where J 1 , . . . , J η are the number of Kautz functions used in each orthonormal projection of the Volterra kernels, the random process (θ, i 1 , . . . , i η ) ∈ Θ×Z η + → B η (θ, i 1 , . . . , i η ) represents the η-order random Volterra kernel, expanded in the orthonormal basis, and (θ, n j ) ∈ Θ × Z + → ψ η,ij (θ, n j ) represents the random version of the i j -th Kautz filter related with the η-order random Volterra kernel. The Kautz functions are modeled as random processes, because their definition 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 written

y(θ, k) ≈ ∞ η=1 J1 i1=1 . . . Jη iη=1 B η (θ, i 1 , . . . , i η ) η j=1 l η,ij (θ, k) , (14)
where the random process (θ, k)

∈ Θ × Z + → l η,ij (θ, k) is a simple filtering of the deterministic input signal u(k) by the random Kautz function ψ η,ij , i. e., l η,ij (θ, k) = V -1 ni=0 ψ η,ij (θ, n i )u(k -n i ) , (15)
for V = max{J 1 , . . . , J η }.

Then, the coefficients of the kernels can be estimated, considering Monte Carlo simulations and the least squares method. Figure 1 shows a flowchart of the Volterra kernels identification considering the random Kautz functions and the Monte Carlo simulations. So, considering each stochastic realization θ, the matrix Γ can be completed with the input signal filtered l η,ij (θ, k) and the vector y with the experimental output signal y(θ, k)

Φ = (Γ T Γ) -1 Γ T y , (16)
where Φ has the terms of the orthonormal kernels B η , in each realization θ. The procedure is repeated until the Monte Carlo convergence is achieved.

The stochastic Kautz parameters

In order to create a robust identification, it is necessary to provide some type of certification, i. e., an envelope of reliability for the nominal values of the model parameters. Such certification can be obtained by a stochastic model, where probability distributions, instead of scalar values, are identified 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

θ ∈ Θ → ζ n (θ) ∈ R, θ ∈ Θ → ω n (θ) ∈ R.
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 θ ∈ Θ → ξ η (θ) ∈ R, θ ∈ Θ → ω η (θ) ∈ R.
So, the random system response is used to define the Kautz poles values. Thus, the Kautz parameters, related with the first Volterra kernel, can be approximated as natural frequencies and damping ratios of the equivalent linear 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 definition of the Kautz parameters 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 define linear approximations between modal parameters and the high order Kautz parameters

ω 1 (θ) = ω n (θ) , ξ 1 (θ) = ζ n (θ) , ω 2 (θ) = p 1 ω n (θ) , ξ 2 (θ) = p 2 ζ n (θ) , (17)
ω 3 (θ) = p 3 ω n (θ) , ξ 3 (θ) = p 4 ζ n (θ) ,
where the relationships p 1 , p 2 , p 3 and p 4 are estimated minimizing the error function min

∆ J = N k=1 [ŷ(k) -y(k)] 2 (18)
where y(k) is the deterministic discrete-time system response, ŷ(k) is the deterministic model response and the vector ∆ contains the coefficients of the relationships between modal and Kautz parameters. Once the coefficients were defined, 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.

Damage detection based on Volterra series

The stochastic Volterra series can be used as a mathematical model to approximate the system response, first in the healthy condition to define 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 version 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 identified 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 identified 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 confidence. In this work, two different approaches are shown, based on the kernels coefficients and on the linear and nonlinear contributions to the total system response, in each structural condition. The approaches are described as follow.

Damage detection based on kernels coefficients

The first approach is based on the coefficients of the Volterra kernels identified. This methodology does not consider the direct influence of the Kautz functions, because only the coefficients expanded in the orthonormal basis are considered. The influence of the random Kautz functions is taken into account only in the process of model identification. Considering Eq. ( 16) and the Volterra series truncated in the third order term, the reduced order kernels coefficients, for each Monte Carlo realization, are allocated as

Φ =   Φ 1 Φ 2 Φ 3   , Φ 1 =          B 1 (1) B 1 (2)
. . .

B 1 (J 1 )          , Φ 2 =                          B 2 (1, 1) B 2 (1, 2) . . . B 2 (1, J 2 ) B 2 (J 2 , 1) . . . B 2 (J 2 , J 2 )                          , Φ 3 =                                    B 3 (1, 1, 1) B 3 (1, 2, 1) . . . B 3 (1, J 3 , 1) B 3 (1, 1, 2) . . . B 3 (1, J 3 , J 3 ) . . . B 3 (J 3 , J 3 , J 3 )                                    , (19)
remembering that J 1 , J 2 and J 3 are the number of Kautz functions used, considering the first, second and third kernels, respectively. It is expected that the main information about the high order kernels is allocated in the main diagonal, so it is considered as monitoring parameters only the coefficients positioned in this region. Then, the coefficients with higher contribution can be used to represent the model identified

λ 1 =          B 1 (1) B 1 (2) . . . B 1 (J 1 )          , λ 2 =          B 2 (1, 1) B 2 (2, 2) . . . B 2 (J 2 , J 2 )          , λ 3 =          B 3 (1, 1, 1) B 3 (2, 2, 2) . . . B 3 (J 3 , J 3 , J 3 )          , λ 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 identified in reference condition

λ 1 =          B 1 (θ, 1) B 1 (θ, 2) . . . B 1 (θ, J 1 )          , λ 2 =          B 2 (θ, 1, 1) B 2 (θ, 2, 2) . . . B 2 (θ, J 2 , J 2 )          , λ 3 =          B 3 (θ, 1, 1, 1) B 3 (θ, 2, 2, 2) . . . B 3 (θ, J 3 , J 3 , J 3 )          , λ nl = λ 2 λ 3 . (21)
Additionally, the index m is used to represent the four different indicators calculated. So, the general notation λ m is used from here, with m = 1, 2, 3, or nl, depending on the number of the kernel considered in the analysis. The damage detection can be summarized as follows 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.

Damage detection based on kernels contribution

With the aim of comparison between the results, a different approach is proposed, based on the contribution of the Volterra kernels identified, to the total response. In this situation, the convolution between the input signal filtered 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

y 1 (k) ≈ J1 i1=1 B 1 (i 1 ) l i1 (k) , (22)
y 2 (k) ≈ J2 i1=1 J2 i2=1 B 2 (i 1 , i 2 ) l i1 (k) l i2 (k) , (23)
y 3 (k) ≈ J3 i1=1 J3 i2=1 J3 i3=1 B 3 (i 1 , i 2 , i 3 ) l i1 (k) l i2 (k) l i3 (k) , (24)
y nl (k) = y 2 (k) + y 3 (k) . ( 25
)
where y 1 (k), y 2 (k), y 3 (k) and y nl (k) are the linear, quadratic, cubic and nonlinear contributions, respectively. Again, 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

y 1 (θ, k) ≈ J1 i1=1 B 1 (θ, i 1 ) l i1 (θ, k) , (26)
y 2 (θ, k) ≈ J2 i1=1 J2 i2=1 B 2 (θ, i 1 , i 2 ) l i1 (θ, k) l i2 (θ, k) , (27)
y 3 (θ, k) ≈ J3 i1=1 J3 i2=1 J3 i3=1 B 3 (θ, i 1 , i 2 , i 3 ) l i1 (θ, k) l i2 (θ, k) l i3 (θ, k) , (28)
y nl (θ, k) = y 2 (θ, k) + y 3 (θ, k) . ( 29
)
As mentioned before, the index m is used to represent the five different indexes calculated. So, the general notation y m is used from here, with m = 1, 2, 3, or nl. The damage detection is very similar to the methodology proposed using kernels coefficients and can be summarized as follows 1. Identification 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 (y m ); 3. Identification of the Volterra model in an unknown condition; 4. Calculation of the kernels contribution to the total model response, in an unknown condition (y m ); 5. Comparison between the contributions estimated in unknown and reference condition, if the new model belongs to the set of reference models, the structure is classified as healthy.

The difference here is that using the kernels contribution, the Kautz functions have an influence on the process, whereas using only the coefficients, they do not. Based on that, the variations in the frequency of oscillation and the damping coefficients 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.

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 (λ m or y m ), the classical Mahalanobis squared distance can be calculated [START_REF] Worden | Damage detection using outlier analysis[END_REF]

D m (θ) = [λ m (θ, i) -µ λm ] T Σ -1 [λ m (θ, i) -µ λm ] (30)
or

D m (θ) = [y m (θ, k) -µ ym ] T Σ -1 [y m (θ, k) -µ ym ] (31)
where D m (θ) is the Mahalanobis squared distance in the reference condition, Σ is the covariance matrix and µ (.) is the mean operator. Note that, in the reference condition, D m (θ) 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 [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF][START_REF] Worden | Experimental validation of a structural health monitoring methodology: Part i. novelty detection on a laboratory structure[END_REF] pDm

(d m ) = 1 N h N i=1 K d m -d m,i h , (32)
where pDm (d m ) is an approximation of p Dm (d m ) that is the true density of D m (θ), N is the total number of realizations of the random parameter, d m,i is the ith realization of the random variable D m (θ), K represents the kernel of the transformation, in this work the Gaussian kernel is used, and h is the smoothing parameter that controls the width of the Gaussian kernel chosen. The main difficulty 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 [START_REF] Bowman | Applied smoothing techniques for data analysis: the kernel approach with S-Plus illustrations[END_REF]. The Kernel Density Estimation is used, because the distribution of D m (θ) is assumed as unknown a priori and the large number of Monte Carlo simulations, provides sufficient samples to estimate its density using this approach.

Then, with the density estimated pDm (d m ) it is possible to establish a threshold value for the distribution, which can be used in the damage detection procedure. The threshold value can be defined as [START_REF] Rébillat | Peaks Over Threshold-based detector design for structural health monitoring: Application to aerospace structures[END_REF] Λ m = {d m such that

+∞ dm p Dm (d m ) dd m = β} , (33)
where Λ m represents 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/confidence that each operational application demands. Now, a new model is identified in an unknown condition (healthy or damaged). The indexes are computed considering the new model (λ unk m or y unk m ) and compared with the stochastic model (set of reference models), through the Mahalanobis squared distance

D unk m = [λ unk m -µ λm ] T Σ -1 [λ unk m -µ λm ] (34)
or

D unk m = [y unk m -µ ym ] T Σ -1 [y unk m -µ ym ] (35)
where D unk m is 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

H 0 : D unk m ≤ Λ m , H 1 : D unk m > Λ m , (36)
where the null hypothesis H 0 represents the healthy condition and H 1 the damaged. The methodology used to detect damage is summarized in the flowchart showed in fig. 2, considering the two proposed approaches.

Modeling a breathing crack in an uncertain nonlinear beam

This section describes the studied nonlinear system that can be approximated by a nonlinear single-degree-of-freedom (SDOF) model, a Duffing oscillator [START_REF] Kovacic | The Duffing equation: nonlinear oscillators and their behaviour[END_REF]. 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.

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 positioned in its free extremity, to generate a nonlinear interaction with a magnet (fig. 3). The objective of the setup used is to emulate a hardening nonlinear 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 displacements and its dynamic behavior can be well approximated, for a limited frequency range, by its first mode shape that can be represented as a Duffing oscillator [START_REF] Kovacic | The Duffing equation: nonlinear oscillators and their behaviour[END_REF]

m ẍ(t) + c ẋ(t) + k 1 x(t) + k 2 x(t) 2 + k 3 x(t) 3 = U (t), ( 37
)
where m is the system equivalent mass in [kg], c is the damping coefficient in

[Ns/m], k 1 is the linear stiffness in [N/m], k 2 is the quadratic stiffness in [N/m 2 ],
k 3 is the cubic stiffness in [N/m 3 ], and U (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), ẋ(t) and ẍ(t). The quadratic stiffness was added to emulate the interaction between shaker and beam, to ensure a more realistic description of the system behavior, whereas the cubic stiffness describes the interaction between beam and magnet. The deterministic values of the Duffing 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, considering experimental data measured. The description of the procedure used to approximate the real system considered by the Duffing 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 effect 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 confirm 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 stiffnesses of the Duffing oscillator motion equation. These results show, clearly, the nonlinear behavior of the system for large displacements (condition reached considering 1 N as excitation amplitude).

Damage simulated -A breathing crack model

The methodology showed in section 2 proposes to detect damages, or structural 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, considering that the system has nonlinear behavior before (Duffing oscillator) and after (breathing crack) the damage occurrence, taking into account the data variation. In other words, the nonlinear characteristic of the damage has different nature of the nonlinear behavior of the system in the reference condition. A 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 [START_REF] Chati | Modal analysis of a cracked beam[END_REF][START_REF] Andreaus | Non-linear dynamics of a cracked cantilever beam under harmonic excitation[END_REF][START_REF] Chatterjee | Structural damage assessment in a cantilever beam with a breathing crack using higher order frequency response functions[END_REF][START_REF] Prawin | Nonlinear structural damage detection based on adaptive Volterra filter model[END_REF]. Combining the reference model (Duffing oscillator) with the cracked beam model, the system behavior can be approximated by

m ẍ(t) + c ẋ(t) + F(x(t), t) + k 2 x(t) 2 + k 3 x(t) 3 = U (t), (38)
where F(x(t), t) is the bilinear force

F(x(t), t) = k 1 x(t), if x(t) ≥ 0 α k 1 x(t), if x(t) < 0
where 0 < α < 1 represents the crack severity. If α = 1 the system is in the reference condition (without crack). The breathing crack phenomenon generates a nonlinear behavior with different nature of the cubic effect of the reference model but similar to the quadratic effect of k 2 . This situation makes the application of damage detection procedures difficult. 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 [START_REF] Peng | Crack detection using nonlinear output frequency response functions[END_REF] [START_REF] Peng | Crack detection using nonlinear output frequency response functions[END_REF], [START_REF] Surace | Detecting nonlinear behaviour using the Volterra series to assess damage in beam-like structures[END_REF] [START_REF] Surace | Detecting nonlinear behaviour using the Volterra series to assess damage in beam-like structures[END_REF] and [START_REF] Prawin | Nonlinear structural damage detection based on adaptive Volterra filter model[END_REF] [START_REF] Prawin | Nonlinear structural damage detection based on adaptive Volterra filter model[END_REF], for example, that have used the Volterra series to detect damage, considering the linear behavior in reference condition and the nonlinear effects as a consequence of the crack occurrence.

In the presence of the breathing crack, the nonlinear behavior does not depend 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 presence 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.

Stochastic version of the nonlinear system

It is known that mechanical systems and structures are subjected to uncertainties that have to be considered in the processes of models identification and damage detection procedures, as it was exemplified in section 1. So, taking the uncertainties into account is essential to make robust predictions [START_REF] Cunha | Quantification of parametric uncertainties induced by irregular soil loading in orchard tower sprayer nonlinear dynamics[END_REF]. 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 difficult the application of damage detection procedures [START_REF] Sohn | Effects of environmental and operational variability on structural health monitoring[END_REF]. The main idea is to simulate these variabilities in order to reproduce, in a more realistic way, the difficulties involved in the damage detection 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 identification and damage detection approach proposed in section 2, the linear stiffness (k 1 ) and damping coefficient (c) of Eq. ( 38) are considered as random variables θ

∈ Θ → k 1 (θ) ∈ R, θ ∈ Θ → c(θ) ∈ R.
The variation of k 1 and c reflects directly in the system dynamics, since these parameters influence the fundamental frequency of oscillation and the damping ratio of the system. So, the Eq. ( 38) can be rewritten as

m ẍ(θ, t) + c(θ) ẋ(θ, t) + F(x(θ, t), t) + k 2 x(θ, t) 2 + k 3 x(θ, t) 3 = U (t), ( 39
)
where the bilinear random force F(x(θ, t), t) can be described as

F(x(θ, t), t) = k 1 (θ) x(θ, t), if x(θ, t) ≥ 0 α k 1 (θ) x(θ, t), if x(θ, t) < 0
and the random processes (θ, t) ∈ Θ × R → x(θ, t), (θ, t) ∈ Θ × R → ẋ(θ, t), and (θ, t) ∈ Θ × R → ẍ(θ, 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 [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Cunha | Modeling and Quantification of Physical Systems Uncertainties in a Probabilistic Framework[END_REF]. Considering that the linear stiffness and the damping coefficient can not be negative, we assumed the interval (0, ∞) as the support of these random variables. Also, it was considered that the expected value of k 1 and c are known real numbers µ k1 and µ c . For technical reasons, see Soize (2017) [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF] for details, we also supposed that the expected value of ln k 1 and ln c are finite. Using these conditions as known information, such as done in Cunha and Sampaio (2015) [START_REF] Cunha | On the nonlinear stochastic dynamics of a continuous system with discrete attached elements[END_REF], the principle of maximum entropy says that the probability density function (PDF) of these random variables is given by [START_REF] Silva | Non-parametric identification of mechanical systems by Kautz filter with multiple poles[END_REF] where Z represents the random parameter (k 1 or c), µ Z is the mean value, δ Z is the dispersion, Γ indicates the gamma function, and 1 (0,∞) denotes the indicator function of the interval (0, ∞). This PDF corresponds to a gamma distribution. The stochastic formulation was made assuming, for the lack of better knowledge, that k 1 and c are independent random variables. To compute the propagation of the uncertainties of the random parameters k 1 and c through the model, the Monte Carlo (MC) method [START_REF] Rubinstein | Simulation and the Monte Carlo Method[END_REF] 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 × 10 3 [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 identification and damage detection procedures.

p Z (z) = 1 (0,∞) 1 µ Z 1 δ 2 Z 1 δ 2 Z 1 Γ(1/δ 2 Z ) z µ Z 1 δ 2 Z -1 exp - z δ 2 Z µ Z
To exemplify the difficulty to detect damage in this situation, fig. 8 shows the 99% confidence bands of the system Frequency Response Function (FRF), considering 3 different 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 situation is not recommended, because only conditions of high severity damage can be detected with probabilistic confidence. 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 described, considering the initially nonlinear behavior and data variation. (a) Low level of input (0.1 N).

Damage detection procedure

This section shows the principal results obtained with the application of the damage detection approach proposed, based on the stochastic Volterra series, 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 samples. The Monte Carlo simulations were performed considering 2048 realizations, to ensure the method convergence. The response of the system described in section 3 was simulated considering different 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 conditions) and the simulated data variation. The simulations regarding the refer-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.

Stochastic Volterra models identification

The identification of the stochastic Volterra model was done considering 2048 realizations, for each structural condition, obtained through the simulations performed. Only the first three kernels were used to identify the system, once its nonlinear characteristic in the healthy and damaged condition can be approximated 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 influence 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 define the Kautz parameters, taking into account the relations described in Eqs. [START_REF] Lim | Reference-free fatigue crack detection using nonlinear ultrasonic modulation under various temperature and loading conditions[END_REF]. Table 2 shows the values of the factors p 1 , p 2 , p 3 and p 4 found for the high order kernels, considering the different 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 effect of the crack. The number of functions were chosen as done in [START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF] [START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF], and defined as J 1 = 2, J 2 = 4 and J 3 = 6. With the Kautz poles defined, the random Kautz functions can be obtained and used in the models' identification process. The identification 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 (first mode region). The kernels identification was carried out in two steps, first, 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) [START_REF] Shiki | On the application of discretetime Volterra series for the damage detection problem in initially nonlinear systems[END_REF]. The models were estimated via MC simulation, in a total of 2048 realizations for each structural condition.

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, defined by

conv(N ) = 1 N N n=1 t f t=t0 ||h(θ n , t)|| 2 dt , (41)
where N is the number of MC simulations, ||.|| denotes the standard Euclidean norm and h(θ n , t) represents the n-th realization of the random first kernel or main diagonal of high order kernels, in the time domain. For further details, the reader is encouraged to see Soize (2005) [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF]. The criterion was applied to the first, second and third kernels. Figure 9 shows the results obtained, considering the Volterra kernels identified in the reference condition. It can be seen that the MC convergence is achieved with the number of samples used.

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 confidence. Thus, the same input used to estimate the Volterra kernels was applied to the model. Figure 10 shows the comparison between the 99% confidence 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.2 0 0.2 0.4 0.6 2 2.2 2.4 2.6 2.8 3 3.2 3.4 Velocity [m/s] Time [s]
(a) Complete Signal (α = 1.00).

- Then, a signal with different 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 the visualization of the harmonic components, with 99% of statistical confidence 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 (fig. 11a) and damaged (fig. 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 amplifies 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. 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 filter and compare the nonlinear contributions to the total system response before and after the damage occurrence. In order to exemplify this information, fig. 12 shows the contributions of the first, second and third kernels with 99% of statistical confidence, 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 influenced by the crack behavior. This result shows that the kernels contributions, or kernels coefficients, can be used in the process of damage detection. The variation of the quadratic contribution is bigger, because of the nonlinear nature of the damage. Finally, in the situation of severe damage, the difference between the response in healthy and damaged condition can be visu-ally observed, but at the beginning of the crack propagation, the differentiation is more difficult, 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.

Damage detection

This section exemplifies the two different metrics proposed to detect damage in nonlinear systems, using the stochastic version of the Volterra series. As mentioned 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 stiffness and damping coefficient. The first 2048 realizations in the reference condition were used as training data and the others as test data.

Use of kernels coefficients

The approach proposed, based on Volterra kernels coefficients, showed in section 2.6.1 was applied, combined with the novelty detection showed in section 2.6.3. The multiple models identified in all structural conditions were considered in the calculation of the indexes. Four different indexes were calculated (λ 1 , λ 2 , λ 3 and λ nl ), to compare the results. Figure 13 shows the evolution of the Mahalanobis squared distance applied to the indexes, with the crack propagation. The boxplot method is used to clarify the visualization of the results. The first 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 coefficients 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 coefficients of the second and third kernels. The dispersion of the values is small, because the uncertainties in the modal parameters have lower influence in the estimation of the higher order Volterra kernels coefficients. The dispersion observed in the quadratic index is related to the noise added because the quadratic component has the level of amplitude close to the noise. With the Mahalanobis squared distance of the indexes calculated in the reference condition, and considering their empirical distribution obtained, it is possible to define threshold values, taking into account a probability of false alarms (β) required in each application. In this work, three different values were considered, β = 0.005, 0.01 and 0.02. The results of the hypothesis test applied are shown in fig. 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 coefficients) 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 identification, but, nevertheless, the nonlinear index was able to make difference 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 confidence. With the aim of study the performance of the approach proposed, the receiver operating characteristics (ROC) curve was computed [START_REF] Farrar | Structural health monitoring: a machine learning perspective[END_REF]. This curve presents the relation between true detections and false alarms for different β values. Figure 15 shows the results obtained for the 4 indexes. The linear index has the worst performance and the cubic index an intermediate performance. 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 influence in the coefficients 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.

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 section 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 response. The difference compared with the approach presented in the last section is the influence of the Kautz functions in the process, that introduces the influence of the natural frequency and damping ratio in the approach, through the Kautz poles dependence. The increase of the crack makes the natural frequency 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 influence using this approach, doing the indexes vary more.

Again, the multiple models identified 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, considering a high level of amplitude (1 N), to maximizes the nonlinear influence in the response. Four different indexes were calculated (y 1 , y 2 , y 3 and y nl ), 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 performance to detect the presence of the crack compared with the use of kernels coefficients. The increase of the indexes is related to the variation of the natural 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 calculated in the reference condition, and considering their empirical distribution obtained, it is possible to define threshold values, taking into account the probability of false alarms (β) required in each application. Again, three different values were considered, β = 0.005, 0.01 and 0.02. The results of the hypothesis test applied are shown in fig. 17, for the three values of β used. The performance 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 [START_REF] Farrar | Structural health monitoring: a machine learning perspective[END_REF]. 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 variation 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 influence 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 first kernel contribution.

Comparison between the use of kernels coefficients and kernels contribution

Comparing the results that have been shown in figs. [START_REF] Worden | Experimental validation of a structural health monitoring methodology: Part i. novelty detection on a laboratory structure[END_REF] and 16, it can be noticed that the indexes variation is bigger using the kernels contribution, because 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 coefficients. On the other hand, the capability of the Kautz functions to describe the variation of the frequency of oscillation related with the crack presence makes the approach more sensitive to the presence of the damage, 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. Additionally, comparing the figs. 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 influence 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, like the breathing crack phenomenon studied, the coefficients of the higher order kernels showed to be sensitive to the presence of the damage. In this situation, the approach based on the kernels coefficients is more effective to make difference between the uncertainties and the damages. Therefore, both approaches have shown interesting performance and different 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 considered, 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 coefficients it is interesting.

Final remarks

The paper has presented an approach to detect damages in initially nonlinear systems, considering data variation caused by the uncertainties. The approach is based on a stochastic version of the Volterra series expanded using random Kautz functions. Two different methodologies were presented, considering the Volterra kernels coefficients and the Volterra kernels contribution. The final decision about the state of the system (healthy or damaged) was performed 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 proposed have presented similar performance. As expected, the linear analysis fails when the damage has low severity, because of the influence 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 behavior is considered, in this work, even in the reference condition, which combined with the presence of uncertainties, makes difficult the damage detection process. But, the nonlinear indexes have shown to be able to differ the initial nonlinear behavior and data uncertainties of the crack behavior, with satisfactory performance. 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.

Figure 1 :

 1 Figure 1: Description of the Volterra kernels identification approach used, based on Monte Carlo simulations.

1 .

 1 Identification 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 coefficients, in the reference condition (λ m ); 3. Identification of a new Volterra model in an unknown condition; 4. Construction of the indexes based on kernels coefficients, in an unknown condition (λ m ); 5. Comparison between the coefficients estimated in unknown and reference conditions, if the new model belongs to the set of reference models, the structure is classified as healthy.

Figure 2 :

 2 Figure 2: Description of the damage detection approach based on stochastic Volterra series.

Figure 3 :

 3 Figure 3: Nonlinear system simulated in the reference condition.

Frequency

  Low input (0.1 N).

Frequency

  High input (1 N).

Figure 4 :

 4 Figure 4: Time-frequency diagram of the response for different levels of input amplitude, considering the reference condition (Healthy).

Figure 5 :

 5 Figure 5: Nonlinear system simulated in the damaged condition.

Frequency

Frequency

Figure 6 :

 6 Figure 6: Time-frequency diagram of the response for different levels of input amplitude, considering the crack severity α = 0.9.

Figure 7 :

 7 Figure 7: Distribution of the random parameters k 1 and c.

  High level of input (1 N).

Figure 8 :

 8 Figure 8: FRFs considering different structural conditions and 2 levels of input amplitude. B represents α = 1 (reference condition), B α = 0.98 (damaged condition) and B α = 0.96 (damaged condition).

6 N

 6 (c) Main diagonal of third kernel.

Figure 9 :

 9 Figure 9: MC convergence test applied to the Volterra kernels identified in the reference condition (α = 1.00).

Figure 10 :

 10 Figure 10: Stochastic Volterra model response in comparison with new simulated data, considering a high level amplitude chirp input (1 N) and two structural conditions. The gray box represents the 99% confidence bands, --represents the mean and -• -the simulated data.

  Damaged condition (α = 0.86).

Figure 11 :

 11 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% confidence bands, --represents the mean and -• -the simulated data.

Figure 12 :

 12 Figure 12: Confidence bands of the kernels contribution to the total system response obtained through the Volterra model with 99% of confidence, considering high level of input (1 N). B represents the reference condition (α = 1.00) and B the damaged condition (α = 0.86).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86 Mahalanobis squared distance Crack severity (α)(a) Linear index (λ 1 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86 Mahalanobis squared distance Crack severity (α) (b) Quadratic index (λ 2 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86 Mahalanobis squared distance Crack severity (α) (c) Cubic index (λ 3 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.86 Mahalanobis squared distance Crack severity (α) (d) Nonlinear index (λ nl ).

Figure 13 :

 13 Figure 13: Boxplot of the Mahalanobis squared distance applied to the four indexes calculated, for different levels of crack severity.

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (a) Linear index (λ 1 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (b) Quadratic index (λ 2 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (c) Cubic index (λ 3 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (d) Nonlinear index (λ nl ).

Figure 14 :

 14 Figure 14: Percentage of damage detection obtained using the kernels coefficients, considering different crack severities and probability of false alarms used. B represents β = 0.02, B represents β = 0.01 and B represents β = 0.005.

Figure 15 :

 15 Figure15: Receiver operating characteristics (ROC) curve, obtained using the kernels coefficients. --represents the linear index, -x -the quadratic index, --the cubic index and --the nonlinear index.

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96 Mahalanobis squared distance Crack severity (α)(a) Linear contribution (y 1 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96 Mahalanobis squared distance Crack severity (α) (b) Quadratic contribution (y 2 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96 Mahalanobis squared distance Crack severity (α) (c) Cubic contribution (y 3 ).

  1.00 0.98 0.96 0.94 0.92 0.90 0.88 0.96 Mahalanobis squared distance Crack severity (α) (d) Nonlinear contribution (y nl = y 2 + y 3 ).

Figure 16 :

 16 Figure 16: Boxplot of Mahalanobis squared distance of the kernels contribution calculated for different levels of crack severity.

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (a) Linear contribution (y 1 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (b) Quadratic contribution (y 2 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (c) Cubic contribution (y 3 ).

  0.98 0.96 0.94 0.92 0.90 0.88 0.86 Percentage of detection [%] Crack severity (α) (d) Nonlinear contribution (y nl = y 2 + y 3 ).

Figure 17 :

 17 Figure 17: Percentage of damage detection obtained using the kernels contribution, considering different crack severities and probability of false alarms used. B represents β = 0.02, B represents β = 0.01 and B represents β = 0.005.

Figure 18 :

 18 Figure 18: Receiver operating characteristics (ROC) curve, obtained using the kernels contribution. --represents the linear index, -x -the quadratic index, --the cubic index and --the nonlinear index.

Table 1 :

 1 Duffing oscillator parameters.

	Parameter Deterministic value
	m [kg]	0.26
	c [Ns/m]	1.36
	k 1 [N/m]	5.49 × 10 3
	k 2 [N/m 2 ]	3.24 × 10 4
	k 3 [N/m 3 ]	4.68 × 10 7

Table 2 :

 2 Relation factors between the Kautz and modal parameters, considering the high order Volterra kernels.

	Crack severity Second kernel Third kernel
	(α)	p 1	p 2	p 3	p 4
	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

Acknowledgments

The authors are thankful for the financial support provided from of São 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.