Yaşar Yanık
email: ysryanik@gmail.com

Samuel Da Silva
email: samuel@dem.feis.unesp.com

Americo Cunha
email: americo@ime.uerj.br

Uncertainty quantification in the comparison of structural criterions of failure

Keywords: failure criterions, stress analysis, uncertainty quantification, parametric probabilistivc approach, Monte Carlo method

Uncertainty Quantification is a discipline which deals with the quantitative description and reduction of uncertainties in real world applications. It tries to find out how likely certain results are in a few parts of the system are not precisely known. Furthermore, failure theory is the investigation of predicting the circumstances under which solid materials under the processing of external loads. The Failure Theory is known as different failure criteria such as Von Mises and Tresca which are the most famous of these for certain materials. In this context, this paper is going to show a comparison between Tresca and Von Mises failure criterions, taking into account the underlying uncertainties in the constitutive equations and stress analysis.

INTRODUCTION

The aim behind the failure criteria is to predict or estimate the yield and failure of machine parts and basic individuals. For the most part the Von Mises and Tresca criteria are demonstrated together with little proposal or discrimination between them in all mechanics of materials books. Furthermore, the primary explanation of the Mises criterion is that it shows a critical value of the distortional energy keeped in the isotropic material while the Tresca criterion is that of a crucial value of the maximum shear stress in the isotropic material [START_REF] Burns | The theory of materials failure, by richard m. christensen: Scope: monograph. level: postgraduate, advanced undergraduates[END_REF]. Additionally, most of the predictions that are essential for decision making in engineering, financial aspects, actuarial sciences are made in view of computer models. These models depend on presumptions that could be not in accordance with reality. Therefore, a model can have uncertainties on its forecasts, because of conceivable wrong presumptions made during its originations [START_REF] Cunha | Modeling and Uncertainty Quantification in the Nonlinear Stochastic Dynamics of Horizontal Drillstrings[END_REF][START_REF] Cunha | Modeling and Quantification of Physical Systems Uncertainties in a Probabilistic Framework[END_REF]. Additional varieties may rise from an diversity of sources including the geometry of the issue, material characteristics, limit conditions, starting conditions, or excitations forced on the framework [START_REF] Wojtkiewicz | Uncertainty quantification in large computational engineering models[END_REF]. The parametric probabilistic approach [START_REF] Soize | Stochastic models of uncertainties in computational mechanics[END_REF] comprises in modeling the uncertain parameters of the computational model by random variables, in order to construt a stochastic model to deal with the underlying variabilities. This kind of approach is highly fited and and extremely effective to consider the uncertainties on the computational model parameters. Many studies have been published in this area, see [START_REF] Soize | Stochastic modeling of uncertainties in computational structural dynamics -recent theoretical advances[END_REF] for details.

Once a computational model is contructed, it is necessary investigate how the uncertainties propagate from model parameters to the response. This can be done via Monte Carlo (MC) method, that is the most widely recognized technique for stochastic calculation, because of its simpleness and great factual outcomes. Nevertheless, its computational expense is to a exceedingly, and, as a rule , inhibitive. Luckily the Monte Carlo calculation is effectively parallelizable, which permits its utilization in simulations where the calculation of an individual realization expensive [START_REF] Cunha | Uncertainty quantification through the monte carlo method in a cloud computing setting[END_REF]. This technique generates several realizations (samples) of the random parameters according to their distributions (stochastic model). Each of these realizations defines a deterministic problem, which is solved (processing) using a deterministic technique, generating an amount of data. Then, all of these data are combined through statistics to access the response of the random system under analysis [START_REF] Cunha | Uncertainty quantification through the monte carlo method in a cloud computing setting[END_REF].

In this paper, the quantitication of security factor uncertainties is addressed through a parametric probabilistic approach, using MC as stochastic solver. The goal is to is to verify the effect of the parametric uncertainties in the safety factor obtained with basis on the failure criteria of Mises and Tresca.

DETERMINISTIC FAILURE CRITERIONS

At this section, it is shown that the state of plane tension of the mechanical system of interest is defined by the stress tensor

σ = σ x τ xy τ xy σ y , (1)
where the parameters σ x and σ y are the normal stresses and τ xy is shear stress respectively, being obtained trhough the solution of the eigenvalue problem

|[σ] -λI|υ = 0, (2)
where (λ, υ) is an eigenpair for [σ]. The obtained eigenvalues are equal to the principal stress values σ 1 and σ 2 used in Von Mises failure criterion,

σ vm = (σ 2 -σ 1 ) 2 + σ 2 1 + σ 2 2 , (3)
and Tresca failure criterion

τ tresca = σ 1 -σ 2 . ( 4
)
Thus, the safety factor for Von Mises failure criterion is defined as

F s(V onM ises) = σ e σ vm , (5)
while the safety factor for Tresca failure criterion is given by

F s(T resca) = σ e σ 2 -σ 1 , (6)
where σ e is the yield stress of the material.

SOTOCHASTIC FAILURE CRITERIONS

In this paper, stochastic version of stress matrix is proposed to identify the comparison of structural criterions of failure. The model parameters subjected to uncertainties are described as random variables. Therefore, the system response also becomes a random variable. Therefore, it is shown that the state of plane tension of the mechanical system of interest is defined by stochastic version of the stress tensor

σ = σ x τ xy τ xy σ y (7)
In this part of section, the value of σ y is given as 2 MPa in the Eq.( 8). After giving this value, random values are used by applying uncertainty quantification method considering τ xy and σ x . Afterward, these values were examined in the non-standard experimental distribution graph considering Maximum Entropy Principle [START_REF] Kesavan | The generalized maximum entropy principle[END_REF] and the σ values were obtained. Furthermore, stochastic version of the eigenvalue equation is described as shown below

|[σ] -λI|υ = 0. ( 8
)
In addition, while obtaining unknown parameters in these equations, 1024 samples, obtained via Monte Carlo method, are shown on the plot. According to Figs.3, convergence plot of Tresca failure criterion begins to converge after 200 number of samples. At a value of 1024 number of samples, the graph is approximately stable. Moreover, stochastic version of the equation of safety factor for Tresca failure criterion is determined as shown below

F s(T resca) = σ e σ 2 -σ 1 . ( 9
)
Stochastic version of stress value for Von Mises will be obtained as shown below

σ vm = (σ 2 -σ 1 ) 2 + σ 2 1 + σ 2 2 . ( 10
)
After obtaining stochastic version of stress value for Von Mises σ vm , stochastic version of safety factor for Von Mises failure criterion is determined as shown below

F s(V onM ises) = σ e σ vm . (11)

Maximum Entropy Principle and Convergence Criteria

According to Maximum Entropy Principle, the stochastic equivalent of the approximation shown in Eqs.( 12), ( 13) and ( 14) are written

+∞ -∞ p x (x)dx -1 = 0, (12) +∞ -∞ xp x (x)dx -µ = 0, ( 13
) +∞ -∞ x 2 p x (x)dx -µ 2 -σ 2 = 0, ( 14
)
where p x (x) is probability density function of the random parameter. When these three equations, ( 12), ( 13) and ( 14) are solving, it is also assumed that support, µ (mean) and σ (variance) of the random variables are known.

p x (x) = e -λ0 e -λ1x-λ2x 2 (15)
where λ 0 , λ 1 and λ 2 are the solution the Eqs.( 12), ( 13) and ( 14). Considering Eqs.( 12), ( 13) and ( 14), the values of σ are given as 280 kPa and 120 kPa separetely for solving these equations in stochastic way for σ x and τ xy random variables. Also, µ (mean) and σ (variance) values of the σ x and τ xy are given as 200 kPa and 280 kPa respectively for solving these equations. After giving these values, the variables λ 0 , λ 1 and λ 2 were obtained. These variables λ 0 , λ 1 and λ 2 equal to -75.162, -0.71425 and 0.00178564 respectively for solving these non-linear equations.

Furthermore, the mean-square convergence analysis [START_REF] Cataldo | Probabilistic modeling of a nonlinear dynamical system used for producing voice[END_REF] with respect to independent realizations

F s(T resca) (θ 1 ),..., F s(T resca) (θ n ) of the random variable F s(T resca) is carried out studying the function n -→ Conv(n) defined by Conv(n) = 1 n n j=1 F s(T resca) (θ j ) 2 (16)
Also, non-standard experimental distribution graphs are described for solving the stress tensor considering σ x normal stress and τ xy shear stress shown in Figure 1. The probability density function that was obtained according to the Von Mises and Tresca criterion using Eq.( 6), is interpreted, the blue and pink colors describe Tresca failure criterion and Von Mises failure criterion, respectively. According to Figs. 3, in the range of 2.0-3.5 F s (Safety Factor), the areas where Von Mises and Tresca failure criterion have a common value range. In addition, the closest match between Von Mises and Tresca failure criterion is in the range of 2.7 F s (Safety Factor). Also, there is no common intersection area between Von Mises and Tresca failure criterion in the range of 1.5-2.0 F s (Safety Factor) and 3.5-4.3 F s (Safety Factor), respectively. It is seen that Von Mises and Tresca failure criterion are not different from each other according to Figs. The convergence plot that was obtained according to Tresca criterion using Eq.( 5), is to be interpreted. In the range of 0-10000 number of samples, Tresca failure criterion is irregularly distant from each other according to Figs. 5. Furthermore, in the range of 10000-250000 number of samples, convergence values between Tresca failure criterion was found to approach each other

The Simple Deflection Problem Considering Uncertainty Quantification

The cantilever beam AB is uniform cross section and carries a load F at it is free and end at the point of A. The equation of the elastic curve and the deflection and slope at A were determined as shown below [START_REF] Beléndez | Large and small deflections of a cantilever beam[END_REF].

d 2 v dx 2 = M (x) EI (17)
Using the free-body diagram of the portion AB of the beam according to Figs. 6, if the equation is determined as shown below

M (x) = F (L -x) (18)
if the equation integrate Eq.( 18) considering variable x, the equation will be obtained as shown below

EI dv dx = F Lx - 1 2 F x 2 + C 1 (19)
We now observe that at the fixed end B we have v(0) = 0 and θ = dv/dx = 0. Substituting these values into Eq.( 19) and solving for C 1 , then C 1 is determined as shown below

C 1 = 0 (20)
Figura 6: Deflection Beam Integrating both members of Eq.( 19), then the equation is determined as shown below

EIv = 1 2 F Lx 2 - 1 6 F x 3 + C 2 (21)
But, at B we have x = L , y = 0. Substituting into Eq.( 21), then the equation is determined as shown below

1 2 F Lx 2 - 1 6 F x 3 + C 2 = 0 (22) C 2 = 0 (23)
Carrying the value of C 2 back into Eq.( 21), then the equation of the deflection function is determined as shown below

v(x) = 1 EI 1 2 F Lx 2 - 1 6 F x 3 (24)
The boundary conditions are that the displacement and slope are both zero at clamped and from which the two constant of integration can be obtained. Then, the slope function will be maximum at x = L, the equation of the slope function is determined as shown below

dv dx = 1 EI F Lx - 1 2 F x 2 (25) dv dx = 1 EI 1 2 F L 2 (26)
Afterward, the Hooke-Lamé's Law in Cartesian Coordinates were examined and the σ x values were obtained. In this part of the work, the Hooke-Lamé's Law in Cartesian Coordinates were examined in a deterministic way considering ε x displacement and other variables were eliminated. After that, the eliminated equation of Hooke-Lamé's Law in Cartesian Coordinates is determined as shown in Eq.( 28)

        σ x σ y σ z τ xy τ xz τ zy         =         λ + 2G λ λ 0 0 0 λ λ + 2G λ 0 0 0 λ λ λ + 2G 0 0 0 0 0 0 G 0 0 0 0 0 0 G 0 0 0 0 0 0 G                 ε x ε y ε z γ xy γ xz γ zy         (27)
σ x = λ + 2G ε x (28)
The slope function at x = L will be equal to displacement ε x .

dv dx = ε x = 1 EI 1 2 F L 2 (29)
In the Hooke-Lamé's Law in Cartesian Coordinates, shear (G) and lambda (λ) modulus were determined [START_REF] Landau | Theory of elasticity[END_REF] as shown in Eqs.( 30) and ( 31) for explain to σ x normal stress

λ = νE (1 + ν)(1 -2ν) (30)
Considering Eqs.( 12), ( 13) and ( 14), the values of F max and F min are given as 350N and 150N separetely for solving these equations in stochastic way. Also, µ (mean) and σ (variance) values of the force are given as 250N and 350 N respectively for solving these equations. After giving these values, the variables λ 0 , λ 1 and λ 2 were obtained in Matlab. These variables λ 0 , λ 1 and λ 2 equal to -93.1334, -0.71428 and 0.0014286 respectively for solving these non-linear equations.

Also, non-standard experimental distribution graphs are described for solving the simple deflection problem shown below In addition, according to Figs. 9, these graph shows that comparison between Von Mises and Tresca failure criterion of probability density functions. In addition, the closest match between Von Mises and Tresca failure criterion is in the range of 1.18 F s (Safety Factor). Also, there is no common intersection area between Von Mises and Tresca failure criterion in the range of 0.8-0.96 F s (Safety Factor) and 1.36-1.93 F s (Safety Factor). As a result, Von Mises and Tresca failure criterion has many common intersection areas in the range of 1.1-1.36 F s (Safety Factor) and it is seen that Von Mises and Tresca failure criterion are not different from each other according to Figs. 9 considering uncertainty quantification method for F

FINAL REMARKS

The using of Von Mises and Tresca failure criterions can help the investigation of predicting the circumstances under which solid materials under the processing of external loads. The results obtained in the present paper have shown the comparison between Von Mises and Tresca failure criterions by taking into the consideration with Monte Carlo Simulations and Maximum Entropy Principle. Based on the results shown, the failure criterions investigated here can be useful in applications for the identification and analysis of stress analysis. Future steps

Figura 1 :Figura

 1 Figura 1: The comparison between generated and theoretical non-standard distrubution

Figura 3 :

 3 Figura 3: Probablility Density function and Histogram for Von Mises and Tresca

3 4. 1

 31 Results and Discussion Considering Two Uncertainties In this part of the section, the random parameters τ xy and σ x are independent in each other. According to the probability density function graph Figs. 4, the F s values were obtained in the Von Mises and Tresca equations, and Von Mises and Tresca failure criterion methods were compared. Furthermore, after obtaining the probability density function graph Figs. 4, a convergence graphic Figs. 5 between Von Mises and Tresca was acquired too. Also, in the range of 2.65 F s (Safety Factor), the plots of probability density function for failure criterions intersect in each other according to the probability density function graph Figs. 4.

Figura 4 :

 4 Figura 4: Probability Density Function and Histogram for Von Mises and Tresca considering two uncertainties

Figura 5 :

 5 Figura 5: Convergence of the Tresca failure criterion considering two uncertainties

Figura 7 :Figura 8 :

 78 Figura 7: The comparison between generated and theoretical non-standard distrubution

Also, moment of inertia of cantilever beam is determined as shown in Eq.( 32) below [START_REF] Landau | Theory of elasticity[END_REF][START_REF] Timoshenko | History of strenght of materials: with a brief account of the history of theory of elasticity and theory of structures[END_REF][START_REF] References Beléndez | Experimental study of the bending of a cantilever beam[END_REF]. In this part of the section, where the parameters of stress matrix are E (elastic module) , I (moment of inertia), F (force) and L (length of beam) respectively according to Eq.( 26) slope function, the value of E and L are given as 210 * 10 9 P a and 1m for solving the equation in a deterministic way

Stochastic Explanation for the Simple Deflection Problem

In this part of the section, non-standard experimental distribution graph is used for solving the stress tensor and the simple deflection problem instead of γ distribution graph. Because, when the negative values are used for normal stress σ x and Force F . Non-standard experimental distribution graph should be used considering Maximum Entropy Principle. After that, random variables τ xy and σ x are considered for stress matrix. The process is repeated until the convergence is achieved. Furthermore, the parameters F for the equation of the slope function are considered as random variables according to simple beam deflection problem. Afterward, these parameters were examined in the non-standard experimental distribution graph considering Maximum Entropy Principle [START_REF] Pavon | On the geometry of maximum entropy problems[END_REF].

Using the free-body diagram of the portion AB of the beam according to Figs. 6, if the stochastic version of equation is determined as shown below

if the equation integrate Eq.( 33) considering variable x, the equation will be obtained as shown below

Integrating both members of Eq.( 32), then the stochastic version of equation is determined as shown below

Carrying the value of C 2 back into Eq.( 35), then the stochastic version of the deflection function is determined as shown below

The slope function will be maximum at x = L, the equation of the stochastic version of slope function is determined as shown below

After that, the eliminated the stochastic version of Hooke-Lamé's equation is determined as shown in Eq.( 39)

The stochastic version of slope function at x = L will be equal to displacement ε x .

ACKNOWLEDGMENTS

The authors would like to thank the financial support provided by CAPES, CNPq, FAPESP, FAPERJ.


Luis G G Villani
email: luisgustavovillani@gmail.com

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

Americo Cunha
email: americo@ime.uerj.br.

Damage detection based on a stochastic version of the Volterra series

This paper aims the damage detection in a clamped-free beam, presenting nonlinear behavior for large displacements, considering the data variation related with the presence of uncertainties. For this, a stochastic version of the Volterra series, expanded through a random version of the Kautz functions, is implemented, using Monte Carlo simulations to estimate the underlying statistics. The measurements were performed in different days and conditions, to consider the data variation and the presence of nonlinear behavior in the analysis. Two different damage indexes are proposed, based on the advantage of the use of Volterra series, related with the capability to separate linear and nonlinear contributions in the total system response. The results have shown that the nonlinear approach is more sensible to the presence of inserted damage (loss of mass), presenting better results than the classical linear approach.

Introduction

Many structures and systems can operate in nonlinear regime of motion, depending on the working conditions that they are subjected. So, they exhibit complex responses with harmonic distortion, jumps, modal interactions, bifurcation, and possible chaos [START_REF] Noel | Nonlinear system identification in structural dynamics: 10 more years of progress[END_REF]. In Structural Health Monitoring (SHM) problems, this nonlinear behavior can be confused with the presence of inserted damage, mainly when the classical approaches, based on linear metrics, are implemented [START_REF] Bornn | Damage detection in initially nonlinear systems[END_REF]. In these situations, more sophisticated methods have to be used, making the approach able to make difference between the intrinsic nonlinear behavior of the system and the presence of damage.

On the other hand, any measurement is uncertain in nature, because of the variability related, for example, with the presence of noise, environment conditions, sensors sensibility, and others [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF]. These uncertainties can reflect in data variation and to influence the model's parameters estimation. Then, this variability can also difficult the application of damage detection procedures, considering real structures and systems. Therefore, it is essential a statistical analysis to ensure the confidence of the obtained results and to reduce the number of false positives in the SHM approaches applied [START_REF] Mao | Statistical modeling of frequency response function estimation for uncertainty quantification[END_REF].

In this sense, this paper proposes the use of a stochastic version of Volterra series as an effective tool to take into account measurement uncertainties in the process of damage detection in a clampled-free beam underlying large displacements, using for this the Monte Carlo method to construct estimations of the statistical quantities of interest. The Volterra series are expanded using random Kautz functions to reduce the number of terms to be identified and the processing time. The damage inserted, associated to change of mass in connections (bolt loosed), is detected based on comparison between linear and nonlinear contributions of the Volterra series in the total response, estimated in reference and damaged conditions. The results have shown that to consider the nonlinear terms in the analysis improves the capability of the approach to detect damage with a lower level of false alarms.

Detecting Damage Using Stochastic Volterra series

In this paper, a stochastic version of the classical Volterra series is proposed as a tool to detect damage in a structure operating in nonlinear regime of motion and that it is subject to data variation. In general, the systems uncertainties can be classified in two main groups, the data uncertainties and the model uncertainties [START_REF] Soize | A comprehensive overview of a non-parametric probabilistic approach of model uncertainties for predictive models in structural dynamics[END_REF]. The model uncertainties are related with the lack of knowledge about the physics of the problem, so, they can be reduced or even eliminated with the increase of knowledge about the system behavior [START_REF] Cunha | Modeling and Quantification of Physical Systems Uncertainties in a Probabilistic Framework[END_REF]. This type of uncertainties is not taken into account in this work. Thus, it is assumed that the Volterra series are able to describe the dynamics behavior of the structure. On the other hand, the data uncertainties are related with the model parameters variation, caused by data variation, environment conditions, noise, and others [START_REF] Soize | Uncertainty Quantification: An Accelerated Course with Advanced Applications in Computational Engineering[END_REF]. In this work, the data uncertainties are considered through the use of a parametric probabilistic approach. In this sense, the model parameters subjected to uncertainties are assumed as random variables or random process, defined on the probabilistic space ( , , ) AEP , where is a sample space, AE is a σ-algebra over , and P is a probability measure. It is assumed that any random variable () YR in this probabilistic setting, with probability distribution () P dy Y on R, admits a probability density function (PDF) () y p dy Y with respect to dy.

Stochastic version of the Volterra series

The discrete-time Volterra series can describe the nonlinear relationship between input and output signals through the generalization of the convolution concept [START_REF] Schetzen | The Volterra and Wiener Theories of Nonlinear Systems[END_REF]. For better understanding about the Volterra series formulation and more details considering the deterministic analysis, some other works are suggested [START_REF] Da | Nonlinear identification in structural dynamics based on Wiener series and Kautz filters[END_REF][START_REF] Sb Shiki | On the application of discrete-time Volterra series for the damage detection problem in initially nonlinear systems[END_REF]. Thus, consider a nonlinear system with a single output, that it is a random process ( , ) ( , ) kk Zy , as a consequence of a single input, that it is deterministic ( )

( ) k u k Z
. The stochastic version of the Volterra series can describe the random output of the nonlinear system as

N N i i nn k n n u k n yH , (1)
where the random process represents the random version of the η-order Volterra kernel and N 1 ,…,N η the number of terms used. In this work the Volterra series will be truncated in the third order term (η = 3), because of the characteristic of the nonlinear system in study, that can be well approximated by a Duffing oscillator [START_REF] Kovacic | The Duffing equation: nonlinear oscillators and their behavior[END_REF]. The main advantage of the use of Volterra series in the description of nonlinear systems is its capability in separate linear and nonlinear contributions in the total response through the kernels convolution

NN nn k n n u k n u k n yH , ( 3
) 3 12 1 2 3 1 11 3 3 1 2 3 1 2 3 0 0 0 ( , ) ( , , , ) ( ) ( ) ( ) N NN n n n k n n n u k n u k n u k n yH .
There are some drawbacks in the use of this formulation, like the high number of terms, the difficulty in the series convergence and the model stability. So, in order to reduce the high number of terms the kernels can be expanded in a orthonormal bases, such as the Kautz functions, that will be considered in this work because of its capability in describe oscillatory systems [START_REF] Wh Kautz | Transient synthesis in the time domain[END_REF][START_REF] Ps Heuberger | Modelling and identification with rational orthogonal basis functions[END_REF]. Then, it can be considered where represents the η-order random Volterra kernel, expanded in the Kautz basis, J 1 ,….,J η are the number of Kautz functions used in each orthonormal projection and the random process represents the random version of the i j -th Kautz filter. It is important to observe that the Kautz functions are assumed as random, because their definition depends of the system dynamics and, in this work, the system output, which is related with its dynamic behavior, is considered as random. Now the random system output can be estimated through multiple convolutions between the random Volterra kernels, represented in the orthonormal basis, and the deterministic input signal filtered by the random Kautz functions where the random process represents the deterministic input signal, filtered by the random Kautz functions, i. e., where, V=max{J 1 ,…,J η }.

( , ,..., ) ... ( , ,..., ) ( , )

j J J ij j ii n n i i n H B É , (4) 1
1 1 1 1 1 1 ( , ) ... ( , ,..., ) ( , ) j J J i j ii k i i k y B l , ( 5
) 1 0 ( , ) ( , ) ( ) jj i V i i i i n k n u k n lÉ , (6) 1 1
( , ,..., ) ( , ,..., )

i i i i ZB ( , ) ( , ) j j i j nn ZÉ ( , ) ( , ) j i kk Zl 4
Then, the estimation of the coefficients can be done considering the least squares method and Monte Carlo simulations. In this regard, considering each realization θ, the input signal filtered ( , )

j i k l
can be used to complete the matrix Γ, and the random experimental output ( , ) k y used to complete the vector Y. Finally, the coefficients estimation is made considering where, Φ is composed by the terms of the orthonormal kernels B η for each realization θ. The procedure is repeated until the Monte Carlo convergence is achieved. Figure 1 shows a flowchart describing the kernels estimation procedure used in this work, considering the Monte Carlo simulations.

Random version of the Kautz functions

The deterministic version of the Kautz functions can be found in [START_REF] Wahlberg | System identification using Kautz models[END_REF]. The functions definition is related with the system dynamics, expressed as natural frequencies and damping ratios. As, in this work, the system response is assumed as random, the modal parameters are also considered as random. So, the Kautz functions are modeled as random process, of generalize can be rewritten as and being the random parameters and defined as

T 1 T () Φ = Γ Γ Γ Y , ( 7
) 22 2 ,2 1 22 [1 ( ) ][1 ( ) ] ( ) ( )[ ( ) 1] 1 ( , ) ( )[ ( ) 1] ( ) ( )[ ( ) 1] ( ) j z z z z z z z bc c b c É b c c b c c , (8) ,2 ,2 1 2 () ( , ) ( , ) 1 ( ) jj z zz b ÉÉ b , (9) () b ()
c
where and represent, respectively, the random Kautz poles in discrete domain, related with each Volterra kernel, and its complex conjugate. The discrete poles can be related with the continues poles where F s is the sampling frequency. Finally, the Kautz poles, for each kernel, are defined as where and are the parameters of the Kautz poles and η represents the number of the Volterra kernel considered in the analysis. When it is cosidered the linear analysis and can be approximated, respectively, by the damping ratio and natural frequency , that are considered random variables as a consequence of the data variation. The definition of the parameters related with the high order Volterra kernels it is more complicated, because they are not exactly the natural frequencies and damping ratios, but combinations of these values. In this work, it is proposed a linear approximation of these parameters, considering the kernels of second and third order where the relationships, p 1 , p 2 , p 3 and p 4 can be determined through the minimization of the error function where is the system response obtained through the Volterra model, is the experimental system response, and Δ the vector with the coefficients p 1 , p 2 , p 3 and p 4 . The optimization is realized only for a representative data, and not in all Monte Carlo simulations, this methodology reduces the computational costs. So, in this work, the system's mean response is used in the optimization procedure. Then, the Kautz functions poles are obtained through the modal parameters estimated in each Monte Carlo simulation and the coefficients estimated.

[ ( ) ( )] () 1 ( ) ( ) WW b WW , ( 10
) () ( ) exp s F S W , (11) 2 ( )
( ) ( ) ( ) 1 ( ) j S Ø ÿ ÿ Ø , ( 12
) 21 ( ) ( ) n p Ø°, 22 ( ) ( ) n p ÿÿ , 33 ( ) ( ) n p Ø°, 34 ( ) ( ) n p ÿÿ . ( 13
) 1 min J [ ( ) ( )] N k y k y k , ( 14
) () W () W () Ø () ÿ 1 () Ø 1 () ÿ () n °() n ÿ ˆ() yk () yk

Damage detection approach based on novelty detection

The damage detection is done based on the variation of the Volterra kernels contributions in the total system response and novelty detection. As mentioned before, the main advantage of the Volterra series is the capability to separate linear and nonlinear contributions through the kernels convolution, and this advantage allows the comparison between the contributions before and after the damage occurrence. So, this characteristic allows the analysis of which term of the response is more sensible to the presence of damage. In this sense, two damage indexes can be calculated, in reference condition, considering the contributions obtained through the Eq. ( 3)

1 ( , ) ( , ) lin ref ref kk yy , (15) 23 ( , ) ( , ) ( ,
)
nlin ref ref ref k k k y y y ,
where is a matrix N s × k representing the linear contribution, composed by only the first kernel contribution, and is also a matrix N s × k representing the nonlinear contribution composed by the sum of the second and third kernels contributions. Additionally, the index m is used to represent the two different indicators calculated. So, the general notation is used from here, with m = lin or nlin, depending on the number of kernels considered in the analysis. N s is the number of Monte Carlo simulations considered and k represents the discrete-time. It is observed that the contributions of the total response obtained, in reference condition, are random processes, estimated with statistic confidence and not simple curves. Therefore, to detect system variations related with damage, is needed to determine if the new contributions obtained through models identified in unknown system condition are contained in these random processes. To create a classification damaged/healthy state a distance-based method to outliers' detection is used [START_REF] Wu | A measure of dna sequence dissimilarity based on mahalanobis distance between frequencies of words[END_REF]. In this sense, two different damage indexes, in reference condition, based on the standardized Euclidian distance between vectors, can be calculated as where and are, respectively, the linear and nonlinear indexes calculated in reference condition, that are random variables, and V is a diagonal matrix N s × k in which the diagonal elements are the standard deviations of the columns of or . The indexes are calculated in a stochastic sense, because the model is stochastic and the model response is assumed as a random process. It is expected that the statistic confidence of the model does the approach able to differ data variation and damage. With the indexes identified, their probability density function (PDF) can be estimated considering Kernel Density Estimation approach [START_REF] Worden | Experimental validation of a structural health monitoring methodology: Part i. novelty detection on a laboratory structure[END_REF][START_REF] Bw Silverman | Density estimation for statistics and data analysis[END_REF] where represents an approximation of the true PDF of , is the ith realization of the random variable , K is the kernel used in the transformation, a

1T 1 ( ) [ ( , ) ( , )] [ ( , ) ( , )] s m N ref ref ref ref ref m m m m n k n k V k n k Õ y y y y , ( 16
) , 1 1 ˆ() s m m i ref m m ref ref N ref i s pK N h h à , (17)

Ã

Gaussian kernel in this work, and h is a smoothing parameter. The main difficulty in Kernel Density Estimation is the choice of optimal value to the smoothing parameter, then 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 index distribution is assumed as unknown a priori and the large number of Monte Carlo simulations provides sufficient samples of the index to estimate its density with this approach.

Then, a threshold of the reference condition can be defined based on the empirical distribution estimated [START_REF] Rébillat | Peaks over threshold-based detector design for structural health monitoring: Application to aerospace structures[END_REF] where β is the probability of false alarms considered and Λ m is the threshold estimated. Now, with the model constructed in the reference condition and the threshold calculated, new Volterra models can be calculated in unknown conditions (healthy or damaged) and compared with the stochastic reference system response where represents the kernels contributions of the Volterra series identified in the unknown condition. So, the hypothesis test can be applied to determine if the system is in healthy or damaged condition where the null hypothesis H 0 represents the healthy condition and H 1 the damaged. Figure 2 shows a flowchart of the damage detection approach proposed in this work, based on a stochastic version of the Volterra series.

Experimental damage detection in a nonlinear beam

The experimental setup considered was the same used in [START_REF] Sb Shiki | On the application of discrete-time Volterra series for the damage detection problem in initially nonlinear systems[END_REF] and [START_REF] Villani | Damage detection in an uncertain nonlinear beam[END_REF], composed by a clamped-free beam with dimensions of 300×18×3[mm 3 ] and a magnet positioned near to the free end, to generate a nonlinear interaction with the beam (Figure 1). To emulate the presence of damage a bolted connection is positioned 150 [mm] from the free end with four nuts with 1 [g] each one, so the reference condition considers the 4 nuts and the damage is simulated through the loss of mass. Table 1 shows the structural conditions considered in the analysis. The signals acquisition was done considering sampling frequency of 1024 Hz and 4096 samples and using a m+p VibPilot data acquisition system. The structure was excited using a MODAL SHOP shaker (Model Number: K2004E01) attached 50 mm from the clamped, considering three different levels of voltage amplitude 0.01 V (low level), 0.10 V (medium level) and 0.15 V (high level). Finally, a vibrometer laser Polytec (Model: OFV-525/ -5000-S) and a Dytran load cell (Model: 1022V) are used to measure the velocity at the free end of the beam and the force excitation, respectively.

{ such that ( )d } ref m m m m ref m ref ref ref m p à , (18) 1T 1 [ ( ) ( , )] [ ( ) ( , )] s m
N unk unk ref unk ref m m m m n y k n k V y k n k yy , (19) 0 1 :
: The stepped sine test was performed to show the nonlinear behavior of the system, close to the first vibration mode of the structure, when the high level of input is applied. Figure 3 shows the results obtained, it is observed the jump phenomena, characteristic of nonlinear systems, when 0.15 V of voltage amplitude is applied, so the structure presents nonlinear behavior for large displacements. Additionally, to consider the data variation in the procedure, the measurements were repeated 160 times in different days during 2 weeks, and then, a white noise with 5% of the signals' RMS value was randomly added to the data, generating 2048 synthetic experimental realizations to be used in the Monte Carlo simulations. To exemplify the difficulty of detect damage in this situation, figure 4 shows the envelope of the FRFs curves, with 99% of confidence, in reference condition, and different levels of damage. It is observed that it is not possible to differentiate the damaged and healthy conditions, mainly when the level of the damage is low (3 nuts). Thus, to use only the inspection of the output in time or frequency domain is not enough to correctly detect the structural state of the system.

unk mm unk mm H H , (20) () unk m yk

Volterra models estimation

The estimation of the Volterra kernels was done in two steps, first a low level (0.01 V) chirp, varying the frequency from 10 to 50 Hz with rate of 10 Hz/s, was applied to estimate the first kernel (linear), and then, a high level (0.15 V) chirp, varying the frequency from 10 to 50 Hz with rate of 10 Hz/s, was applied to estimate the second and third kernels (nonlinear). The procedure is repeated 2048 times, to achieve the Monte Carlo convergence and to construct the stochastic reference model. The reference model was validated in time and frequency domains, considering as input the same chirp signal used in the identification process and a single tone sine with frequency close to the first vibration mode of the structure (~23 Hz) and high level of amplitude (0.15 V). Figures 6 and7 show the results obtained in the validation process. It is observed that the probability bands of the model, with 99% of confidence, can describe the nonlinear system behaviour, in all frequencies components, and the data variation. The main advantage of Volterra series approach is the capability to separate linear and nonlinear contributions through the kernels convolution. This characteristic can be used to detect damage in the system. So, to exemplify, figures 8 and 9 shows the kernels contributions in the total response, considering the responses showed in figures 6 and 7.

The first and third kernels have high contribution in the total response in time domain, the second kernel has low contribution because the output is approximated symmetric.

In frequency domain, the linear kernel contributes only to the first frequency component, the second kernel contributes only to the second harmonic and the third kernel contributes to the first and third harmonics. This occurs because of the mathematical definition of the third kernel that it has dependence of the first kernel.

Damage detection

After the identification of the set of reference models, the stochastic model obtained is used in the damage detection process. Then, new Volterra models are identified in each structural condition, and again, to consider the data variation, 2048 realizations were considered. So, the methodology explained in section 2.3 is applied, considering the two different indexes proposed (linear and nonlinear). The results were obtained considering a chirp input signal, with high level of voltage amplitude and varying the frequency from 10 to 50 Hz with rate of 10 Hz/s. Figure 10 shows the probability to detect damage using the proposed indexes, considering the 5 different structural conditions and 3 values of β. The nonlinear index has better performance, because it is more sensible to the presence of the damage. The linear index cannot classify all damage conditions as damaged. Finally, the receiver operating characteristics (ROC) curve was computed (Farrar 2012). This curve shows the relation between false alarms and true detections obtained considering different values of β and structural conditions. Figure 11 shows the results obtained to the two indexes proposed. The nonlinear index has better performance than the linear index, with lower number of false alarms and higher number of true detections. Therefore, although the damage emulated (loss of mass) it has linear characteristic, the nonlinear metric proposed present better results, because the nonlinear terms of the response are more sensible to the presence of damage, without high influence of the data variation.

Conclusions

This paper has proposed an approach based on a stochastic version of the Volterra series, estimated based on Monte Carlo simulations, to detect damage in a nonlinear system considering data variation related with uncertainties. The results have shown that the nonlinear index proposed has better performance than the classical linear approach, despite the linear characteristic of the damage emulated (loss of mass). This occurs, probably, because the third kernel's mathematical definition is related with the linear kernel, how can be observed on the cubic contribution in the first harmonic of the system's output PSD (Fig. 9), doing the cubic contribution more sensible to the presence of damage. These results prove the advantage in using nonlinear metrics when the monitored structure presents nonlinear behavior. Besides that, the problem is more difficult to solve when uncertainty in the date is assumed and it can change or fluctuate by chance the metrics. In order to separate these effects, the proposed method in this paper has shown adequate to allow identifying the source of nonlinearities or presence of damage in an uncertainty scenario.

Figure 1 .

 1 Figure 1. Flowchart of the random kernels coefficients estimation based on Monte Carlo simulations.

Figure 2 .

 2 Figure 2. Flowchart of the damage detection approach based on Volterra kernels contributions and novelty detection.

Figure 4 .

 4 Figure 4. Stepped sine test considering different levels of input amplitude.

Figure 5 .

 5 Figure 5. FRF representing the data and damage variations.

Figure 6 .

 6 Figure 6. System output, considering a chirp input with high level of amplitude (0.15 V).

Figure 7 .

 7 Figure 7. PSD of system output, considering a sine input with high level of amplitude (0.15 V).

Figure 8 .

 8 Figure 8. Kernels contributions in the total output, considering a chirp input with high level of amplitude (0.15 V).

Figure 9 .

 9 Figure 9. PSD of kernels contributions, considering a sine input with high level of amplitude (0.15 V).

Figure 10 .

 10 Probability of damage detection for different structural conditions.

Figure 11 .

 11 Receiver operating characteristics (ROC) curve.

Table 1 . Structural conditions.

 1

	State	Ref	Dam I	Dam II	Dam III	Rep
	Structural condition	4 nuts (reference)	3 nuts	2 nuts	1 nut	4 nuts (repair)

Acknowledgements

The authors are thankful for the financial support provided from of São Paulo Research Foundation (FAPESP), grant number 2012/09135-3, 2014/02971-6 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.


Christelle Eyraud
email: christelle.eyraud@fresnel.fr

Rodolphe Vaillon

Amélie Litman

Jean-Michel Geffrin

Olivier Merchiers

Polarization effects in 3D vectorial induced current reconstructions

In tomography algorithms, the complex amplitude scattering matrix corresponds to the input parameter. When considering three-dimensional targets, the scattering matrix is now containing vectorial information. Thus, this scattering matrix might be calculated with various polarization projections. Moreover, when dealing with experimental data, we are almost every time faced with truncated data. We focus here on the impact of selecting parts of the amplitude scattering matrix elements versus others and in particular on the influence of the polarization choices on the imaging results. In order to better apprehend the physical content associated to each polarization term, the study is conducted with a simple vectorial induced current reconstruction algorithm allowing to reconstructed qualitative maps of the scene. This algorithm is applied on scaled models of aggregates combined with experimental scattered fields acquired in the microwave frequency range.

Introduction

Electromagnetic wave probing is an interesting tool to obtain the physical features of unknown targets (position, shape, size, complex permittivity). Indeed, these characteristics can be retrieved from the measurements of its scattered field by an inverse procedure. In practice, there exists some difficulties in particular when dealing with measurements. Indeed, as the scattered field can not be measured on an entire surface enclosing the target -due to mechanical constraints or/and due to the measurement time limitation, the data are truncated [START_REF] Bucci | Electromagnetic inverse scattering: retrievable information and measurement strategies[END_REF]. Moreover, all the polarization cases of the scattered field are rarely measured in a multistatic configuration. The measurements are also ineluctably disturbed by noise, which does not a priori follow an uniform or a white Gaussian distribution [START_REF] Geffrin | Optimization of a Bistatic Microwave Scattering Measurement Setup: From High to Low Scattering Targets[END_REF]. It seems thus interesting to understand and analyze which measurements are the most relevant ones. The final goal is to reduce the number of experimental data points while guaranteeing the best reconstructions.

In this contribution, we focus our efforts on the polarization aspect. The polarization of the electromagnetic wave contains useful information as studies in radar polarimetry [START_REF] Cloude | A review of target decomposition theorems in radar polarimetry[END_REF], [START_REF] Roueff | Cramer-Rao lower bound analysis of vegetation height estimation with random volume over ground model and polarimetric SAR interferometry[END_REF], [START_REF] Bellez | A Full Polarimetric Bistatic Radar Imaging Experiments on Sets of Dielectric Cylinders Above a Conductive Circular Plate[END_REF] or light scattering communities [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF], [START_REF] Merchiers | Microwave measurements of the full amplitude scattering matrix of a complex aggregate: a database for the assessment of light scattering codes[END_REF] and references therein [START_REF] Videen | Photopolarimetry in Remote Sensing[END_REF] point out. In inversion procedures, the different polarizations are generally considered to contain the same amount of information and hence, are treated as a whole. In fact, reconstruction results seem to be greatly influenced by the choice of the field components as shown by the results obtained in [START_REF] Eyraud | 3D-aggregate quantitative imaging: experimental results and polariza-tion effects[END_REF]. The present study is performed by calculating the amplitude scattering matrix using two classical polarization projections. The influence on the reconstruction results of each polarization projection will therefore be investigated. Moreover, in order to better apprehend the physical content associated to each term of the amplitude scattering matrix, we have decided to employ a simple induced current reconstruction algorithm which is rendered vectorial as the polarization states of the transmitter and the receiver are taken into account. The reconstructions are made on scaled models of aggregates [START_REF] Vaillon | A new implementation of a microwave analog to light scattering measuremen device[END_REF] with experimental scattered fields acquired in the microwave frequency range [START_REF] Merchiers | Microwave analog to light scattering measurements on a fully characterized complex aggregate[END_REF]. These targets fulfill neither the Born or the Rytov approximation and thus only qualitative maps of the scene will be provided.

The paper is organized as follows. The scattered field is expressed in Part 2. Part 3 is devoted to the analysis of the available spectral information. In Part 4, amplitude scattering matrices are calculated using two conventions, the spherical convention often used in the microwave community and the Bohren and Huffman convention [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] often used in the light scattering community. The induced current reconstruction algorithm, which uses these amplitude scattering matrices as input is described in Part 5. The reconstruction results from the measurements are provided in Part 6 as well as the discussion on the impact of the polarization projection on the reconstructed target. Concluding remarks follow in Part 7.

Far-field scattered field

The field scattered by a target, enclosed in a domain Ω and illuminated with a monochromatic electromagnetic wave, is calculated in a receiving zone Γ, at each receiver position r ∈ Γ, by means of the observation equation which can be written as

E s (r) = Ω G(r, r )χ(r )E(r )dr (1)
A free space configuration is considered, thus G is the dyadic free-space Green function between the object zone Ω and the receiver zone Γ. E and χ (χ(r ) = k(r ) 2 -k 2 0 ) are respectively the field and the contrast in the Ω domain (with k 0 the wavenumber in vacuum and k(r ) the wavenumber in zone Ω).

When the far-field conditions are realized [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF], that is,

k 0 r 1, r r , k 0 r 2r 1, ∀r ∈ Γ, ∀r ∈ Ω(2)
The dyadic free-space Green function in spherical coordinates (Figure 1(a)) can also be approximated by:

G(r, r ) ≈ e jks.(r-r ) 4πr [I -e r ⊗ e r ] ( 3
)
where k s is the scattered wavevector, I the identity operator, a ⊗ b the tensor product between the vector a and the vector b and (e r , e θ , e φ ) are the unit vectors of the spherical basis.

Ewald sphere

If the transmitter and the receiver can be placed all around the Ω zone, on a sphere Γ of radius r, all the information in the far-field is reached for a considered frequency and the Ewald sphere is totally covered [START_REF] Kak | Principles of Computerized Tomographic Imaging[END_REF]. This Ewald sphere corresponds to the location, in the spectral domain, of the vector K = k s -k i , where k i (resp. k s ) corresponds to the incident (resp. scattered) wavevector. Unfortunately this configuration is difficult to obtain in practical situations. Indeed, all geometrical positions may not be allowed due to various considerations (measurement time, mechanical constraints, ...). Figure 1 describes the geometrical configuration associated to the experimental setup of the C.C.R.M. -Marseille employed in the present work for acquiring the data sets (see section 6.A). The origin O of the spherical coordinate system (r, θ, φ) is located at the center of the target zone Ω [START_REF] Zwillinger | CRC Standard Mathematical Tables and Formulae[END_REF]. The transmitter is placed on the sphere surrounding the target, its position is described by (r T , θ T , φ T ) with θ T ∈ [0, 2π] and φ T ∈ [φ T b , φ T u ] (where φ T b and φ T u are the limit φ T -angles). The receiver, placed on the same sphere (r R = r T ), has a position described by (r

R , θ R , φ R ) with θ R ∈ [θ R b , θ R u ] (where θ R b
and θ R u are the limit θ R -angles) and with φ R = 90 • . The associated Ewald sphere which describes the spectral information contained by the scattered field is represented in Figure 2. From a global point of view, this sphere is relatively well-filled. The lack of the backscattering zone -due to the limits of the θ R angle -is balanced The plane (k x 0k y ) is relatively wellfilled. Nevertheless it can be noticed that it contains a larger amount of low spatial frequency values than of high frequency ones. In the (k x 0k z ) and (k y 0k z ) planes, parts of the disks are missing, i.e., the lower and the upper parts along the k z -axis are truncated and only low spatial frequency values are present. This is due to the constrained measurement geometry, i.e., the angles between 0 and φ T b (30 • in the case of Figure 2) are not accessible by the transmitter in the positive part of the z-axis and the angles between φ T u and π (50 • in the case of the Figure 2) are excluded in the negative part. This explains also the asymmetry with respect to the k z = 0 axis. In the following, the discussion on the polarization impact is restricted to this specific Ewald sphere topology.

Amplitude scattering matrix

In the far-field, for a non-absorbing medium, the complex amplitude scattering matrix S, as defined in [START_REF] Mishchenko | Scattering, Absorption, and Emission of Light by Small Particles[END_REF], describes the realtion between the components of the incident field vector and the components of the scattered field vector. With this definition, the scattered field vector E s is linked to the incident field vector E i as:

E s = S E i (4)
This scattering matrix is determined by the incidence and the scattering directions as well as by the properties of the target, but it also depends on the choice of the polarization basis. Two main conventions are usually used, the first one -called here the spherical convention -is well-spread in the microwave community while the second one -here called the Bohren and Huffman convention -is favored by the optics community.

4.A. Spherical convention

For a given transmitter/receiver couple, the scattering matrix can be constructed with a spherical convention, using as basis the (u T θ , u T φ ) vectors for the incident wave and the (u R θ , u R φ ) vectors for the scattered wave [START_REF] Michshenko | Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation[END_REF], [START_REF] Merchiers | Microwave analog to light scattering measurements on a fully characterized complex aggregate[END_REF] (Figures 1(b) and 3(a)). As the far-field conditions are assumed, the field vectors have no radial components and are then totally determined by their components along u R θ and u R φ . Accordingly, the element S p,q of the scattering matrix corresponds to a scattered field measured along the u R q vector for an incident field E i polarized along the u T p vector (for p, q ∈ {φ, θ}) and the scattering matrix is written as: This convention is particularly convenient when the fields are measured with several incident and scattering directions.

E s θ E s φ = S θθ S θφ S φθ S φφ E i θ E i φ (5)

4.B. Bohren and Huffman convention

In the previous definition, the elements are not linked to the scattering plane, i.e, the plane containing both the k i incident wave vector and the k s scattered wave vector (see Figure 3(b)). For a given transmitter/receiver couple, the complex amplitude scattering matrix can be constructed using another basis related to this specific scattering plane, as proposed in the Bohren and Huffman (BH) definition [START_REF] Bohren | Absorption and Scattering of Light by Small Particles[END_REF] which is commonly used in the light scattering community. The scattered field components are then linked to the incident ones in the far-field region as follows:

E s E s ⊥ = S S ⊥ S ⊥ S ⊥⊥ E i E i ⊥ (6)
where E (resp. E ⊥ ) corresponds to the field component parallel (resp. perpendicular) to the scattering plane and the element S p,q of the scattering matrix corresponds to a scattered field measured along the u R q vector for an incident field E i polarized along the u T p vector (for p, q ∈ {⊥, }).

Induced current reconstructions based on scattering matrices

In a classical first-order diffraction tomography algorithm, the scattered field vector and the incident field are along the same direction and the Born approximation or the Rytov approximation is used in Equation (1). This equation becomes a scalar equation and the contrast map is then retrieved via an inverse Fourier transform (see for example [START_REF] Sung | Optical diffraction tomography for high resolution live cell imaging[END_REF], [START_REF] Gurg | Relation between computed tomography and diffraction tomography[END_REF], [START_REF] Devaney | Mathematical foundations of imaging, tomography and wavefield inversion[END_REF]).

Here, a different approach is considered: (i) the entire vectorial nature of the fields is taken into account, and (ii) no linearized assumption is made. We only assume that the induced current vector J(r ) = χ(r ) E(r ) is collinear with the polarization state of the incident field u T p (for p ∈ {φ, θ, ⊥, }). The scattered field expression becomes:

E s (r) ≈ e jks.r 4πr [e θ ⊗ e θ + e φ ⊗ e φ ]u T p Ω e -j(ks-k i ).r J p (r )e -jk i .r dr (7)
The vector u T p describing the polarization state of the incident field is not constant for all the illumination directions. For a given choice of polarization (p ∈ {θ, φ, ⊥ , }), it can be written in the form

u T p = f r p e r + f θ p e θ + f φ p e φ ( 8
)
where the terms f r,θ,φ p are functions of the θ T , φ T , θ R and φ R angles. Due to the specific geometry of the configuration under investigation, the scattered field polarization vector is generally not collinear with the incident field one. The vector u R q , which represents the projection on the selected polarization (q ∈ {θ, φ, ⊥, }) given by the receiver is:

u R q = g θ q e θ + g φ q e φ ( 9
)
where the terms g θ,φ q are functions of the θ T , φ T , θ R and φ R angles. For a selected polarization case (p, q), the element of the scattering matrix S p,q (where p ∈ {θ, φ, ⊥, } and q ∈ {θ, φ, ⊥, }) corresponds then to the scalar product between the vectors u T p and the vector u R q :

S p,q = E s (r).u R q (10)

≈ e jks.r 4πr u T p .u R
q Ω e -j(ks-k i ).r J p (r )e -jk i .r dr

The magnitude of the induced current |J p,q | created by a target at the position r -in a lossless embedding medium -is then linked to the measurement of the element S p,q of the scattering matrix by the following equation :

|J p,q (r )| ≈ A | S p,q (K)| u T p .u R q , if u T p .u R q = 0 (11)
where 1/A = | e jks.r 4πr | is a constant term and S p,q (K) is the 3D inverse Fourier transform of S p,q (r ).

In this work, we only use measurement points with collinear incident and scattered polarization vectors. The S θφ , S ⊥ and S ⊥ elements are thus not exploited herein as they correspond to situations where the incident and the scattered polarization vectors are orthogonal. Nevertheless, these elements are non null only when the target is changing the polarization of the impinging wave. Finally, the number of elements in the complex amplitude scattering matrix which will be considered in the imaging algorithm may vary according to the polarization state. Table 1 summarizes the percent of measurements used in each case.

Reconstruction results from measurements

In this part, the induced current reconstruction algorithm is assessed using scattering matrices determined from experimental data.

6.A. Experimental configuration

To get all possible information on an object in the farfield, the field should be measured on the entire sphere enclosing the target. The spherical experimental setup in the anechoic chamber of the C.C.R.M. -Marseille allows to obtain a configuration similar to the one described in Figure 1. This setup has been used by the researchers of the Institut Fresnel to perform the measurements.

6.A.1. Geometrical arrangement

For a target at the origin of the coordinate system of the spherical setup (Figure 4), the scattered field can be measured on a portion of the sphere enclosing this object. Indeed, the transmitting antenna can move on a circular arch, the target can rotate around the z-axis -which allows to create a situation which is similar to a fixed target and a source displaced over several vertical arches around the target. The receiver is also moved in the horizontal plane with a circular movement centered on the target. This setup was already described in several publications (see for example [START_REF] Geffrin | Continuing with the Fresnel database: Experimental setup and improvements in 3D scattering measurements[END_REF]). As the distance between the center of the setup and the transmitting antenna or the receiving antenna is 102 λ, we are in agreement with the far-field conditions (equation ( 2)) as long as the radius of the minimum sphere enclosing the target has a diameter smaller than 10 λ.

6.A.2. Experimental parameters

The geometrical parameters are chosen in agreement with the maximum spatial bandwidth of the scattered field to measure the field in a non-redundant way [START_REF] Bucci | Electromagnetic inverse scattering: retrievable information and measurement strategies[END_REF]. Indeed, 99 positions of the source are used in this study. Each source has a position (r T , θ T , φ T ) where

φ T ∈ [φ T b = 30 • : 10 : φ T u = 130 • ] and θ T ∈ [θ T b = 20 • : 40 : θ T u = 340 • ]. The field is measured at 53 re- ceiver positions (r R , θ R , φ R = 90 • ) where θ R ∈ [θ R b = θ T + 50 • : 5 : θ R u = θ T + 310 • ].
The measurements were performed for a wavelength comprised between λ = 15 mm and λ = 18.7 mm, but all the results here are for λ = 16.7 mm. Measurements are calibrated to refer to the scattered field associated to an incident field having a unit magnitude and a null phase at the center of the setup [START_REF] Van Den Berg | Blind shape reconstruction from experimental data[END_REF], [START_REF] Eyraud | Validation of a 3D bistatic microwave scattering measurement setup[END_REF]. Data processing is required to extract the maximum of information [START_REF] Geffrin | Optimization of a Bistatic Microwave Scattering Measurement Setup: From High to Low Scattering Targets[END_REF] from the measurements. In particular all the fields are processed to avoid drift phenomena [START_REF] Eyraud | Drift correction for scattering measurements[END_REF].

6.A.3. Polarization

As the two antennas are linearly polarized pyramidal horns (ARAMWH-1826B), the target is illuminated with a linearly polarized wave and one single polarization component of the scattered field is measured at each measurement point. The different polarization cases are obtained by combining two different emitting antennas plus a rotation of the receiving antenna. The experimental data are thus obtained by measuring for each incident field parallel to u T φ (and then u T θ ), the scattered field vectors which are parallel to u R φ and to u R θ (see Figure 1 for notations).

6.B. Targets

We have considered two targets which are clusters of spheres. The spheres of the first object have diameters which are comparable to the wavelength and in theory it should be possible to separate them with our algorithm. On the contrary, the spheres of the second object are much smaller and can not easily be separated with such an algorithm. The first target is a pyramid of four 1.5 λ in diameter spheres of PMMA (ε r = 2.6) (Figure 5(a)). Its total size is equal to 3 λ. It has a density of approximatively 50% as compared to its minimum covering sphere. The sec- . This object was made by stereo lithography, which is a sequential deposit of layers of liquid acrylate resin (ε r = 2.85 + j0.06) [START_REF] Vaillon | A new implementation of a microwave analog to light scattering measuremen device[END_REF]. The diameter of each sphere is equal to 0.3 λ. Its total size is equal to 5 λ. It has a density of approximatevely 3% as compared to its minimum covering sphere. Due to their size and their permittivity contrast with respect to the background, the two targets do not fullfil the validity condition of the Born or the Rytov approximation [START_REF] Wang | A revisit to the validity of Born approximation in high frequency scattering problems[END_REF].

6.C. Reconstruction parameters 6.C.1. Investigation zone

The dimension of the investigation zone is directly linked to the discretisation step in the spectral domain. As the scattered field is not measured with an uniform step along each component of K, the field has to be interpolated before performing the inverse Fourier transform. For the present measurement configuration, the spectral domain step is set equal to 2 leading to an associated investigation zone corresponding to a cube with a side length of 30 λ. As the voxel width in the reconstructed map is directly linked to the maximum values in the spectral domain, the voxel width is bigger along the z-axis (0.65 λ) than along the y-axis (0.27 λ) and along the x-axis (0.28 λ).

6.C.2. Normalized maps

Each reconstructed map contains qualitative information -position and shape -on the target. We only consider magnitude maps and normalize them to obtain values between 0 and 1. A normalized map M n is thus obtained from an initial map M with

M n = |M | -min(|M |) max(|M |) -min(|M |) (12)

6.C.3. Estimate assessment

To have a quantitative assessment of the reconstruction, we use three criteria. Before calculating these criteria, it is necessary to determined the real position of the target mass center d * = (x d , y d , z d ) to avoid bias. Using the reconstructed maps, we determined the real position of this target mass center using a 3D cross correlation between the reconstructed map and the expected one as explained in [START_REF] Eyraud | 3D-aggregate quantitative imaging: experimental results and polariza-tion effects[END_REF]. The first criterion is a comparison between the true and reconstructed normalized maps using a classical L 2 norm.

C L2 = ||M true n -M reco n || L2 (13)
We have shown in [START_REF] Eyraud | 3D-aggregate quantitative imaging: experimental results and polariza-tion effects[END_REF] that the usually adopted criteria L 1 and L 2 -norms are not very well suited to evaluate the accuracy of the reconstruction of complex targets as their values are particularly dependent on the reconstruction of the empty voxels. We thus use two other criteria, i.e., the correlation C c and the cross correlation C cc coefficients between the two maps:

C c = Cov(M reco n , M true n ) V ar(M reco n ) V ar(M true n ) ( 14
)
C cc = r (M true n (r + d * ))(M reco n (d * )) r (M true n (r + 0))(M true n (0)) (15)
where Cov(u, v) is the covariance between u and v, V ar(u) the variance of u, r the position of the considered cell and d * the real position of the target determined as explained in the following. The C cc coefficient quantifies the similarity between the two maps by only comparing the voxels containing some material. The more similar the two maps are, the closer to 0 the criteria C L2 is and the closer to 1 the criteria C cc and C c are. These criteria are calculated considering a reduced (6 × 6 × 6)λ 3 cubic investigation zone centered at 0. Indeed, due to the relative number of voxels without material in the calculation, the considered zone dimension affects the C L2 criterion.

6.D. Reconstructions in the different polarization cases

In this part, results are presented and then discussed.

6.D.1. Pyramid of spheres

The true position of the center of gravity of the target was found to be at d * = (0.12, 0.12, 0.68) λ. This position is very close to the expected one (0, 0, 0.68) λ.

The reconstructions using the amplitude scattering matrix and the spherical convention are plotted in Figures 6, 8 and 7. In each polarization case, the presence of the target can be seen and its size can be determined on these maps. Nevertheless, the quality of the imaging results is strongly different from a chosen polarization case to another. The φφ and φθ cases allow to obtain rather good results, while the θθ case is clearly less accurate. The kind of information contained in each reconstruction is also different according to the polarization cases. Indeed with the S φφ element, the boundaries of the spheres are better reconstructed than the internal parts. On the contrary, with the S φθ element, the internal parts of the spheres are better reconstructed than the boundaries. Looking at the criteria values, the L 2 criterion gives very little information as this criterion is the same (equal to 0.002) for all polarization cases. The values of the two other criteria are given in Table 2. The C c and C cc criteria show that the best cases are the φφ and φθ cases, the latter being slightly better due to the good reconstruction of the internal part of the spheres.

Let us now consider the reconstruction obtained with the scattering matrix elements under the Bohren and Huffman definition (Figures 9 and10). The reconstruction with the S ⊥⊥ element is very close to the one obtained with the S φφ element. This is confirmed by the C cc and C c criteria (Table 2). The result with the S element is not very good and is similar to the one obtained with the S θθ element.

6.D.2. Aggregate of spheres

The real position of the aggregate mass center using the correlation criterion was found at (0, 0.18, 0.31) λ and the expected position was (0, 0, 0.41) λ. This target is a "flat" target, i.e, it has a small thickness along the zaxis as it can be seen in Figure 11. The reconstruction results are plotted at the altitude z = 0 for the different polarization cases in Figures 1213and the criteria values are given in Table 3. The localization and the object global dimensions are well determined for all the polarization cases. As expected, none of the cases allows to distinguish the spheres due to the use of a linear imaging algorithm. In the same way as for the previous object, the reconstruction results differ from a polariza- tion case to another. Indeed, the reconstruction in the θθ case is the worst according to the C c and C cc criteria. The imaging results obtained with the S φφ and S ⊥⊥ elements are very similar -it is also confirmed by the criteria values -and all the "branches" of the target are reconstructed in these cases. The S element provides a better reconstruction than the S θθ element, where all the "branches" of the target are reconstructed, but the boundaries of the reconstruction are blurred. In the S φθ reconstruction result, some "branches" are missing and this leads to lower criteria values, but the others "branches" are correctly reconstructed.

6.D.3. Discussion

For the two targets, the reconstructions from the S φφ and the S φθ elements are rather close to the real object and the reconstructions from the S θθ element is clearly less accurate. With the Bohren and Huffman polarization definition, the S ⊥⊥ element leads to better results than the S element. The results from the S ⊥⊥ case are similar to the ones derived from the S φφ case. Nevertheless, it can be noticed that the S ⊥⊥ element is the only one taking into account all the measurements and this can be of interest depending on the target under test.

The different parts of the target are not reconstructed in a similar way from one polarization case to another even if the polarization of the incident wave is the same. This is clearly visible while comparing the reconstructions from the S φφ and the S φθ cases for the two targets. To better understand this behavior, we have plot-ted the magnitude of the mean of the induced currents calculated in the target zone at the altitude z = 0.45 λ (Figure 14) -corresponding to the middle of the three spheres of lower altitude -for the pyramid of spheres. These currents were obtained using a volume integral algorithm [START_REF] Eyraud | Microwave imaging from experimental data within a Bayesian framework with realistic random noise[END_REF], for source positions between φ T = 30 • and φ T = 120 • with a 30 • step, considering all the θ T angles, and for a single transmitter polarization u T φ . We can see that the magnitudes of the components along e x and e y have stronger values in the central part of the spheres than on the boundaries, in opposition to the components along e z which have an inverse behavior. Thus, the projection of the resulting scattered field along u R φ or u R θ will provide a different informative content. These maps consolidate the results shown in Figures 6 and8, where different parts of the target were reconstructed according to the chosen polarization reception orientation.

Conclusion

In this paper, we have considered the amplitude scattering matrix calculated with two classical polarization projections. The same imaging procedure, a vectorial induced current reconstruction algorithm, has then been used, while the input components varied according to the projections choices.

Our imaging process was tested using the amplitude scattering matrix calculated from experimental fields measured in an anechoic chamber for two clusters of spheres, a pyramid of 4 spheres with a radius of 0.75 λ and an aggregate of 74 spheres with a radius of 0.15 λ. The results show that the polarization choice has a strong influence on the imaging result. With a spherical configuration, the S ⊥⊥ element using the Bohren and Huffman convention or the S φφ and S φθ elements using the spherical convention lead to the best results. A suitable choice of polarization will therefore provide different kinds of information contents, in particular the capability to focus on either the inner or the outer details of the target substructure.

In our linear algorithm, far-field conditions are required and the induced current vector is assumed to be collinear with the incident field polarization vector. Under these assumptions, three-dimensional qualitative maps of the scene are retrieved in real-time with a lowmemory requirement. Depending on the studied target, in particular if the target can create strong depolarization effects, it might be of interest to pursue the investigation with a non-linear imaging algorithm.

Fig. 1 .

 1 Fig. 1. Geometrical configuration.

Fig. 2 .

 2 Fig. 2. Ewald sphere using the specified configuration with φ T b = 30 • , φ T u = 130 • , θ R b = θ T + 50 • and θ R u = θ T + 310 • (each point corresponds to a transmitter-receiver couple): (a) 3D visualization and (b), (c) and (d) cross-sections

Fig. 3 .

 3 Fig. 3. For a given couple of transmitter and receiver positions (source : (r T = 102 λ, θ T = 100 • , φ T = 60 • ), receiver (r R = 102 λ, θ R = 150 • , φ R = 90 • ), representation of the basis vectors for the source (red) and for the receiver (green) (the k i and ks vectors are plotted in black (--)) under (a) the spherical convention, (b) the Bohren and Huffman convention.

Fig. 4 .

 4 Fig. 4. Experimental setup: photography of the anechoic chamber.

Fig. 5 .

 5 Fig. 5. Dielectric targets: Pyramid of 4 spheres (sphere diameter: 1.5 λ) and (b) Aggregate of 74 spheres (sphere diameter: 0.3 λ)

S

  φφ S θθ S φθ S ⊥⊥ S C cc 0.47 0.27 0.53 0.46 0.28 C c 0.73 0.68 0.71 0.72 0.72

Fig. 6 .

 6 Fig.6. Maps reconstructed at several altitudes using the pyramid of spheres scattered fields in the S φφ polarization. The real sphere boundaries are plotted in black.

Fig. 7 .

 7 Fig. 7. Maps reconstructed at several altitudes using the pyramid of spheres scattered fields in the S θθ polarization. The real spheres boundaries are plotted in black.

Fig. 8 .

 8 Fig. 8. Maps reconstructed at several altitudes using the pyramid of spheres scattered fields in the S φθ polarization. The real spheres boundaries are plotted in black.

Fig. 9 .

 9 Fig. 9. Maps reconstructed at several altitudes using the pyramid of spheres scattered fields in the S ⊥⊥ element. The real spheres boundaries are plotted in black.

Fig. 10 .

 10 Fig. 10. Maps reconstructed at several altitudes using the pyramid of spheres scattered fields in the S element. The real spheres boundaries are plotted in black.

Fig. 11 .

 11 Fig.11. 3D visualization (in blue) of the reconstructed maps using the S ⊥⊥ case -the iso-surface threshold value is set to 0.2. The real 74 sphere aggregate is overlaid (in green).

Fig. 12 .

 12 Fig.12. Maps reconstructed at z = 0 λ using different elements of the scattered fields matrix with the spherical convention. The real aggregate boundaries are plotted in green.

Fig. 13 .

 13 Fig. 13. Maps reconstructed at z = 0 λ using different elements of the scattered fields matrix with the Bohren and Huffman convention. The real aggregate boundaries are plotted in green.

  (a)| Jn,x | (b)| Jn,y | (c)| Jn,z | 14. Maps of the normalized induced current magnitude in the target zone at altitude z = 0.45 λ for the pyramid of spheres for a transmitter polarization choice along e φ . Jn,i corresponds to the component along the ei vector (i ∈ {x, y, z}).

Table 1 .

 1 Measurements used with the inversion procedure for the configuration with φT b = 30 • , φ T u = 130 • , θ R b = 50 • and θ R u = 310 • Polarization S φφ S θθ S φθ S ⊥⊥ S

	%	84 55 33 100 46
	Number	4407 2886 1732 5247 2414

Table 2 .

 2 Criteria values for the pyramid of spheres.

Table 3 .

 3 Criteria values for the aggregate.

	S φφ S θθ S φθ S ⊥⊥ S
	L2 0.009 0.008 0.014 0.008 0.009
	C cc 0.70 0.58 0.62 0.70 0.67
	C c 0.51 0.47 0.40 0.51 0.50