Pierre Jehel

Pierre Léger

Adnan Ibrahimbegovic

Structural seismic fragility analysis of RC frame with a new family of Rayleigh damping models

Keywords: Damping, Inelastic time-history analysis, Structural fragility, Fiber beam element, Reinforced concrete frame structure, Earthquake

Structural seismic vulnerability assessment is one of the key steps in a seismic risk management process. Structural vulnerability can be assessed using the concept of fragility. Structural fragility is the probability for a structure to sustain a given damage level for a given input ground motion intensity, which is represented by so-called fragility curves or surfaces. In this work, we consider a moment-resisting reinforced concrete frame structure in the area of the Cascadia subduction zone, that is in the South-West of Canada and the North-West of the USA. According to shaking table tests, we first validate the capability of an inelastic fiber beam/column element, using a recently developed concrete constitutive law, for representing the seismic behavior of the tested frame coupled to either a commonly used Rayleigh damping model or a proposed new model. Then, for each of these two damping models, we proceed to a structural fragility analysis and investigate the amount of uncertainty to be induced by damping models.

Introduction

Decision makers are interested in seismic risk analyses for predicting the postearthquake situation in a given geographical region, so as to anticipate the human, social and economical impact of a major earthquake. For building structures, seismic risk assessment requires three main steps. i) A seismic hazard analysis has first to be performed. It can be either deterministic or probabilistic. In the latter case, the seismic hazard is often expressed as an intensity measure -often the peak ground acceleration (P GA) -with a certain probability of being exceeded in a certain number of years. ii) Then, the seismic fragility of the building considered has to be estimated: it corresponds to the conditional probability P ij = P [DI ≥ DI i | IM = IM j ] of the building to sustain a given damage index DI i for a unique -or a set of -given input ground motion intensity measures IM j . The probability to attain a damage index DI i can then be computed as

P i = j P [DI ≥ DI i | IM = IM j ] • P [IM = IM j ]. (1)
These damage indices have to be related to building performance requirements. iii) Finally, the exposure of the buildings and populations has to be determined.

Structural seismic fragility analysis thus is a key step in the overall earthquake risk management process. This task is commonly achieved by constructing fragility curves from inelastic time-history analyses that take into account the variability in the seismic input motion alone or in both the input motion and the structural model. The main ingredients for fragility analyses are: i) A set of seismic timehistory records representative of the seismic hazard in the geographical region of interest for the project; ii) An inelastic structural model along with a damping model; iii) A mapping between damage indices and structural performance levels; and iv) Statistical tools to analyze fragility curves.

In [START_REF] Howard | Generation of hazard-consistent fragility curves[END_REF], Hwang and Huo present a methodology for constructing fragility curves accounting for uncertainties in the seismic, site, and structural parameters. On the one hand, 8 scenario earthquakes corresponding to different P GA, annual exceedance probabilities, magnitudes, and source-to-site distances are considered; for each scenario, 50 samples of ground motion time histories are generated using a seismological model that takes into account uncertainties in seismic and soils parameters. On the other hand, 50 samples of each of 6 random structural parameters are generated and then combined using the Latin Hypercube sampling technique to eventually generate 50 samples of the inelastic structural model. Then, for each earthquake scenario, the 50 ground motion samples are combined with the 50 structural samples to establish 50 earthquake-site-structure samples. Finally, for each earthquake scenario, 50 values of the damage index are computed from inelastic analyses and fragility curves are constructed.

The preceding approach defines a fully probabilistic approach in the sense that it takes into account uncertainty sources in both the seismic input motions and the inelastic structural behavior. It can however also be worth considering only uncertainties in ground motions. To that purpose, there exists, following Jalayer and Beck [START_REF] Jalayer | Effects of two alternative representations of ground-motion uncertainty on probabilistic seismic demand assessment of structures[END_REF], an alternative to the IM-based approach considered in this work: the probabilistic ground motion time history approach. This latter method is based on a stochastic ground motion model pertaining to seismic source parameters, which has to circumvent the most difficult drawback to be tackled in the IM-based approach, namely whether the selected set of IM thoroughly represent the input signal characteristics. On another hand, the recent work of Rosić et al. [START_REF] B V Rosic | Formulation and computational application of inelastic media with uncertain parameters[END_REF] considers the uncertain structural response of inelastic media and deterministic loading to provide maps of the probability for a component of the stress tensor to reach a given value at a given point.

The concept of fragility curve reduces the vulnerability analysis to the consideration of a unique intensity measure. This limitation has been pointed out in research work where the concept of fragility surfaces emerged. In a study dedicated to the analyze of the limitations of commonly used intensity measures for fragility analysis of single-degree-of-freedom linear and nonlinear systems [START_REF] Kafali | Seismic fragility analysis: Application to simple linear and nonlinear systems[END_REF], Kafali and Grigoriu propose to construct fragility surfaces instead of curves for assessing the seismic performance of nonlinear systems. For a given state of damage in the structure, the proposed fragility surface is the graphical representation of the relationship between the failure probability and the set of intensity parameters (m, r) constituted of the moment magnitude m and source-to-site distance r. This concept of fragility surface is also used in [START_REF] Seyedi | Development of seismic fragility surfaces for reinforced concrete buildings by means of nonlinear time-history analysis[END_REF] by Seyedi et al. who extend it to other intensity measures. Indeed, they construct fragility surfaces that provide the probability for an inelastic reinforced concrete structure to sustain a given inter-story drift ratio, according to the spectral displacement at both eigenperiod T 1 and T 2 . They finally conclude that when dealing with uncertainties propagation, fragility surfaces allow for estimating the variability of structural fragility due to a second IM, which should lead to more accurate seismic risk analyses.

Fragility curves have been used as a baseline to deal with a wide range of issues pertaining to earthquake engineering. In [START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF], Sáez et al. study the effect of considering inelastic dynamic soil-structure interaction on the seismic vulnerability analysis. They construct fragility curves from a very large number of artificially generated input earthquakes. They also introduce the so-called Fisher informa-tion concept which allows for measuring the amount of information contained in the seismic ground motions and thus provide a tool for the statistical analysis of fragility curves. In [START_REF] Saxena | Development of fragility curves for multi-span reinforced concrete bridges[END_REF], Saxena et al. address the issue of assuming identical support ground motion in the analysis of the seismic response of long, multi-span, reinforced concrete bridges. Analyzing fragility curves, they show that considering spatial variation of earthquake ground motions is of first importance. In [START_REF] Popescu | 3d effects in seismic liquefaction of stochastically variable soil deposits[END_REF], Popescu et al. construct fragility curves to present the results of deterministic and both 2D and 3D stochastic analyses of the seismic liquefaction potential of saturated soil deposits. In the context of design rules assessments, Lagaros [START_REF] Lagaros | Probabilistic fragility analysis: A tool for assessing design rules of rc buildings[END_REF] computes fragility curves to analyze the seismic performance of multi-story RC buildings designed according to Greek and European building codes. The fragility curves are plotted from 10,000 simulations based on Monte Carlo techniques to take into account both uncertainty in the seismic signal and in key structural parameters for assessing structural stability.

There is a likely source of uncertainty in inelastic seismic time history analyses which is only rarely considered in fragility analyses, namely damping added to the inelastic structural model so as to introduce in the simulations an amount of energy dissipation coming from inelastic mechanisms that are not explicitly accounted for in the structural model. In [START_REF] Howard | Generation of hazard-consistent fragility curves[END_REF], the critical viscous damping ratio added for the seismic analyses of a reinforced concrete frame building is described by the uniform distribution restricted to the range 2% -4%. In [START_REF] Bruce | Earthquake risk assessment of building structures[END_REF], the structural damping is described by a lognormal distribution with a mean of 2% and a coefficient of variation of 0.62 for the purpose of steel frame fragility analysis. In both these works, Rayleigh damping is added and the uncertainty pertaining to additional damping thus stems from the critical damping ratio, not from the damping model type -Rayleigh in this case. Because it has been shown that it can be difficult to control the amount of damping generated by common Rayleigh damping models throughout inelastic time history analyses [START_REF] Léger | Seismic-Energy Dissipation in MDOF Structures[END_REF][START_REF] Hall | Problems encountered from the use (or misuse) of Rayleigh damping[END_REF][START_REF] Finley | Unintended consequences of modeling damping in structures[END_REF], the main purpose of the work presented in this chapter is to investigate the likely amount of uncertainty introduced by the damping model, in the context of fragility analysis.

To that aim, we proceed as follows. In the next section, we first present a reinforced concrete moment-resisting frame structure -simply referred to as "RC frame" throughout this chapter -which was tested on a shaking table. We use the corresponding experimental data as a reference to validate the developments that we present all along this chapter. Then, we detail in section 3 the numerical model that we use to perform inelastic seismic time-history analyses of the RC frame. The issue of modeling damping is discussed and a new family of Rayleigh damping models is proposed. Results from simulations performed with both a "classical" Rayleigh damping model and the proposed new damping model are compared to the shaking table test results so as to validate the capability of the proposed combination of hysteretic with additional damping models for representing the behavior of the RC frame. The proposed new family of Rayleigh damping models can rely on a physical background which often lacks to commonly used Rayleigh damping models. In section 4, we proceed to the selection of a set of real seismic ground motion records compatible with the seismic hazard in the Cascadia subduction zone. A seismic fragility analysis of the RC frame in this geographical region is then carried out in section 5: fragility curves, along with their statistical analysis, as detailed in [START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF], are constructed for every damping models considered so as to investigate the amount of uncertainty these latter could bring in structural fragility analyses.

RC frame tested on a shaking table

The test structure considered throughout this chapter is represented in Fig. 1. It was designed at a reduced scale of 1/2 according to the provisions of the National Building Code of Canada [START_REF]National building code of canada[END_REF] and of the Canadian concrete standard [START_REF]Design of concrete structures for buildings[END_REF]. The structure was assumed to have a nominal ductility, which corresponds to a force reduction factor R = 2 to compute the design base shear. The various assumptions and parameters used in the design of the this structure can be found in [START_REF] Filiatrault | Seismic performance of ductile and nominally ductile reinforced concrete moment resisting frames. i. experimental study[END_REF]. Four inverted U-shape concrete blocks attached in each span of the beams were used to simulate concentrated gravity loads from framing joints. The centers of gravity of the added masses were computed such that they coincide with the center of gravity of the beams. Service cracks were induced by these added masses. The total weight of the frame was 95kN. The fundamental period T 1 of the structure with added masses was measured at 0.36s in a free-vibration test. Mode 1 excites 91% of the total mass of the structure and when mode 2 is also considered approximately all the mass is excited.

The structure was assumed to be located in a seismic zone 4 in Canada, as depicted in the 1985 seismic zoning map of the National Building Code of Canada. The seismic hazard in this zone is such that peak horizontal ground acceleration between 0.16g and 0.23g is likely to be observed with 10% probability of exceedance in 50 years. Such seismic zones can be found in Western, Eastern and Northern Canada. The ground motion record that was selected for the test program corresponds to the N04W component of the accelerogram recorded in Olympia, Washington on April 13, 1949. Fig. 2 presents the feedback record measured during the test initially scaled to a peak ground acceleration P GA = 0.21g as well as the corresponding elastic response spectrum with 5% viscous damping ratio.

Seismic inelastic time history analysis of the RC frame

The set of equations of motion for the discretized structure is written as:

M d(t) + C(t) ḋ(t) + F R (t) = -M∆ü g (t) (2)
where d(t) is the vector containing the nodal displacements, M is the mass matrix, C(t) is the damping matrix, F R (t) is the inelastic resisting forces vector, Ɵ g (t) is the vector of the rigid body acceleration induced by the ground displacement u g (t).

In this section, we first present the inelastic structural model used to compute F R (t), then we define two damping models: a "classical" Rayleigh model C 1 (t) and a new model C 2 (t), and we finally validate the capability of the two models for representing the behavior of the RC frame presented in the previous section.

Fiber beam/column element

The inelastic structural model is based on a fiber frame element suitably implemented in the framework of a displacement-based formulation so that it can integrate the uniaxial concrete behavior law recently developed by the authors [START_REF] Jehel | Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications[END_REF] and briefly presented in this section. This constitutive model is capable of Figure 2: Acceleration time history recorded on the shaking table during the test and corresponding elastic response spectrum with a critical viscous damping ratio of 5%.

representing the main energy dissipative phenomena likely to occur in concrete: appearance of permanent deformation, strain hardening and softening, stiffness degradation, local hysteresis loops, appearance of cracks. Its theoretical development and numerical implementation are based on thermodynamics with internal variables [START_REF] Germain | Continuum thermodynamics[END_REF][START_REF] Maugin | The Thermodynamics of Nonlinear Irreversible Behaviors -an Introduction[END_REF] and on the finite element method with embedded strong discontinuities [START_REF] Garikipati | A study of strain localization in a multiple scale framework -the one-dimensional problem[END_REF][START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF][START_REF] Oliver | Theoretical and computational issues in modelling material failure in strong discontinuity scenarios[END_REF]. FEAP [START_REF] Taylor | FEAP: A Finite Element Analysis Program, User manual and Programmer manual[END_REF] is the finite element program used for the numerical implementation of the developments presented in this section.

Enhanced kinematics:

The first ingredient of this model is the definition of an enhanced kinematics that takes strong discontinuities into account. This is done, as depicted in Fig. 3 by writing the displacement field u(x, t) as the sum of a continuous displacement ū(x, t) -that is the displacement as it would be in the absence of strong discontinuity -and of displacement jumps ūi (t) in sections Γ i of the solid domain Ω :

u(x, t) = ū(x, t) + n dis i=1 ūi (t)H Γ i (x) (3)
where

H Γ i (x)
is the Heaviside's function which, for a left-to-right oriented domain, is null on the left side of the discontinuity Γ i and unity on its right side.

With the hypothesis of small transformation, we have the following expression Figure 3: Construction of an enhanced displacement field u(x, t) as the sum of a continuous displacement ū(x, t) and of a displacement jump ūi (t)H Γ i (x) pertaining to discontinuity Γ i .

for the normal strain field:

ǫ xx (ū, ūi , t) = ∂ ū(x, t) ∂x + n dis i=1 ūi (t)δ Γ i (x) (4)
where δ Γ i (x) is the Dirac's function.

Stored energy function:

When the Lagrange's variational principle is for instance chosen to derive the governing equations of the system (see Eq. ( 2) for the corresponding discretized form), the internal potential energy U int as to be written:

U int (ū, α, t) = Ω ψ (ū, α, t) dΩ = n c f ib f =1 Ω c f ψ c (ū, α c , t) dΩ c f + n s f ib f =1 Ω s f ψ s (ū, α s , t) dΩ s f ( 5
)
where n c,s f ib is the total number of concrete or steel fibers, Ω c,s f is the volume of the fiber, ψ c,s is the stored energy function for concrete or steel which depends on the continuous displacement field ū(x, t) and on the set of internal variables α c,s . Normal stresses are computed from these functions as

σ xx = ∂ψ ∂ǫ xx . ( 6
)
Set of internal variables:

The set of internal variables α is defined to characterize the evolution of the main energy dissipative -inelastic -mechanisms which develop in the system. These internal variables are the memory of the system. The physical interpretation of each of them is provided in Fig. 4. Note that the constitutive law used here can handle different behavior in compression and tension, and can also reproduce a visco-elastic response (see [START_REF] Jehel | Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications[END_REF] for a full description). Viscosity is not considered in this work and, for the sake of conciseness, we only focus on the compressive part of the behavior law in Fig. 4. The local admissible state of the system is expressed according to criteria functions in the stress-like domain of the set of variables dual to α. When irreversible mechanisms are activated in the structure, internal variables have to be updated and their evolution is governed by the principle of maximum dissipation. From the computational point of view, because we only consider linear hardening and softening laws, there is no need for local iteration when internal variables are updated, except for transitions between hardening and softening regimes, which leads to an efficient resolution procedure.

Inelastic structural model

The finite element mesh and the uniaxial constitutive laws for steel and both confined and unconfined concrete fibers used for the inelastic structural model of the RC frame are shown in Fig. 5. Material behavior laws have been identified to fit experimental monotonic (σ xx -ǫ xx ) curves. The structure is assumed fixed at its base. Rigid end zones are defined to model the beam-to-column connections and rebar slip in surrounding concrete is not represented. These later hypotheses are questionable because the connections exhibit inelastic behavior during the test. This obvious limitations of the structural model has to be reminded when damping model is added in the simulation.

The loading time history consists in two successive phases: i) Static dead load is first applied step by step and then kept constant; ii) The seismic loading is applied. A first validation check of the inelastic structural model is carried out by simulating a free vibration test. Before dead load is applied, the elastic fundamental period of the structure is computed as T ela 1 = 0.28s; then, due to the inelastic behavior of the structure, the elongated period is evaluated as T ini 1 = 0.36s when dead load is completely applied. Both T ela

Proposition of a new family of Rayleigh damping models

Basic definitions:

In the context of inelastic time history numerical analysis, the definition of damping might differ according to the reference cited. On the one hand, in [START_REF]Quantification of building seismic performance factors[END_REF], damping consists in both i) inherent damping resulting from the dissipation of energy by inelastic structural elements, and ii) additional viscous damping added in the simulation to take into account inherent energy dissipation sources not otherwise explicitly considered in the inelastic structural model. On the other hand, in [START_REF]Modeling and acceptance criteria for seismic design and analysis of tall buildings[END_REF], damping is defined as "the portion of energy dissipation that is not captured in the hysteretic response of components that have been included in the model", and it is then suggested in [START_REF]Modeling and acceptance criteria for seismic design and analysis of tall buildings[END_REF] to use "un-modeled energy dissipation" as a more appropriate terminology for damping.

Because, in experimental investigations, measured damping results from all the energy dissipative phenomena, we herein decide to define damping as the combination of both hysteretic damping due to the energy dissipated by all the inelastic phenomena explicitly accounted for in the structural model and additional viscous damping that should be consistent with the inelastic structural model namely, that does not introduce energy dissipation already accounted for in the inelastic structural model.

Problems encountered with Rayleigh damping:

Controlling the amount of additional viscous damping energy dissipated in inelastic time history analyses is a very challenging task [START_REF] Léger | Seismic-Energy Dissipation in MDOF Structures[END_REF][START_REF] Hall | Problems encountered from the use (or misuse) of Rayleigh damping[END_REF][START_REF] Finley | Unintended consequences of modeling damping in structures[END_REF]. This is especially the case for commonly used Rayleigh proportional damping models, that is when the damping matrix is computed, in its most general form, as

C(t) = α(t)M + β(t)K(t), (7)
where K(t) is the tangent stiffness matrix. Several researchers have provided insight in the comprehension of Rayleigh damping regarding the inelastic structural model it is coupled to, have highlighted limitations, and have eventually provided recommendations to cope with them [START_REF] Léger | Seismic-Energy Dissipation in MDOF Structures[END_REF][START_REF] Hall | Problems encountered from the use (or misuse) of Rayleigh damping[END_REF][START_REF] Finley | Unintended consequences of modeling damping in structures[END_REF]. Nevertheless, adding damping and controlling its consistency with the inelastic structural model still remains an issue to be addressed.

Three common phases in seismic response:

We now discuss in a qualitative way the notion of consistency for additional viscous damping. To that purpose, we start by stating that seismic structural response is composed by three main consecutive phases, as illustrated in Fig. 6. Both inelastic structural model and additional damping model must then be capable of representing the salient phenomena corresponding to each of these three phases. Foremost has to be properly modeled what we call here the "key window", namely the time interval within which the major inelastic modifications for structural performance assessment develop. For instance, key mechanisms that control near-collapse structural behavior are listed in [START_REF] Krawinkler | Importance of good nonlinear analysis[END_REF]: degradation of strength and stiffness, and structure P-delta effects. From experimental results, we know that strain rate is another major issue.

A consistent additional damping model should be adapted to each of these three phases as follows: • Phase 1: None or only few incursions in the inelastic domain occur. Energy dissipation in phase 1 thus comes from the friction in the cracks that appeared when applying dead load and from many other mechanisms always present in mechanical systems. When used, visco-elasticity and constitutive laws with local hysteresis [START_REF] Ragueneau | Damage model for concrete-like materials coupling cracking and friction, contribution towards structural damping: first uniaxial applications[END_REF][START_REF] Jehel | Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications[END_REF] in the structural model could account for these energy dissipation sources but a small amount of additional damping usually has to be added. For the RC frame considered in this work, a free-vibration test was carried out after dead load had been applied and a first modal damping ratio of 3.3% was measured [START_REF] Filiatrault | Seismic performance of ductile and nominally ductile reinforced concrete moment resisting frames. i. experimental study[END_REF].

• Phase 2: As the ground motion becomes stronger (at around 8s in Fig. 6), an important amount of seismic energy is imparted to the structure and some parts of the structure then exhibit inelastic behavior. Inelastic struc-tural models are designed to explicitly model part of the numerous inherent nonlinear energy dissipative mechanisms involved in the structural response.

The energy dissipation due to the mechanisms not explicitly accounted for in the inelastic structural model has to be introduced with the additional damping model.

• Phase 3: The structure has suffered irreversible degradations that modified its dynamic properties. Thus, even if the seismic demand is again as low as in phase 1, the energy dissipative mechanisms are different because of frictions in the cracks that appeared within phase 2 or at degraded bound between steel and concrete. Here again, visco-elasticity and behavior laws with local hysteresis [START_REF] Ragueneau | Damage model for concrete-like materials coupling cracking and friction, contribution towards structural damping: first uniaxial applications[END_REF][START_REF] Jehel | Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications[END_REF] in the structural model could account for these damping sources, but it generally has to be completed by additional damping.

Proposition of a new family of Rayleigh damping models:

In the following, two damping models will be used:

• A commonly used Rayleigh model based on tangent stiffness matrix and with two constant coefficients

C 1 (t) = αM + βK(t); (8)
• We propose a new family of models that is directly dependent on both the two key notions in the definition of the three phases introduced above: the capacity of the inelastic structural model to absorb energy and the seismic demand. The model is based on Rayleigh damping with tangent stiffness matrix and with coefficients adapted to each of the three phases:

C 2 (t) = α(t)M + β(t)K(t) (9)
The idea of adapting Rayleigh damping to the capabilities of the inelastic structural model for dissipating energy is present in the use of the tangent stiffness rather than the initial one: it is expected that the choice of tangent stiffness dependent damping will have the main advantage of providing the significant additional source of damping only in the domains/modes that are not accounted for by inelastic model. Such a choice allows to provide the physically based damping phenomena interpretation, which leads to damping coefficients that are easier to identify. The same idea has been further exploited in [START_REF] Tinawi | Seismic safety of gravity dams: from shake table experiments to numerical analyses[END_REF] where 1% viscous damping is added to an inelastic dam model before cracking and 10% after cracking to represent localized high dissipation by friction between crack lips. Another instance is the work presented in [START_REF] Hieu | Structural dynamics of slender ductile reinforced concrete shear walls[END_REF] where viscous damping is added only in the shear wall zones which remain elastic, while no damping is added in the inelastic zones where the structural model is let alone to dissipate the seismic energy.

In spite of its stronger physical background, implementing damping model C 2 (t) is not as straightforward as damping model C 1 (t). First, three sets of Rayleigh coefficients (α p , β p ) p=1,2,3 corresponding to each of the three phases p have to be identified to define the appropriate damping ratios. Second, the instants which delimit the three phases have to be determined, which can be automatically accomplished by the computer program that is capable of detecting the activation of significant inelastic behavior.

Calibration and validation of the models

Seismic inelastic time history analyses of the RC frame have been carried out with the inelastic structural model presented above coupled to either additional damping model C 1 (t) or C 2 (t). The implicit Newmark integration scheme with parameters β = 0.25 and γ = 0.5 is used with a time step of 0.005s. Fig. 7 shows a comparison between the simulated top-displacement and energies time histories and the respective experimental results reproduced from [START_REF] Filiatrault | Seismic performance of ductile and nominally ductile reinforced concrete moment resisting frames. ii. analytical study[END_REF]. Good agreement between simulated and experimental data can be observed. Moreover, there is very good agreement between the hysteretic plus elastic (E H ) and damping (E D ) energy quantities computed with the models proposed here and an analogous Perform3D [START_REF]Perform3D User's manual[END_REF] simulation we carried out for comparison purpose, namely E H ≈ 550N.m and E D ≈ 2250N.m, corresponding to approximately 20% and 80% of the total work done by the structure during seismic motions. Last, the fundamental period of the RC frame in post-earthquake conditions is estimated by a free-vibration test performed at the end of the seismic signal and it is observed that simulated value is T sim 1 = 0.45s whereas the experimental measure comes to T exp 1 = 0.55s. This shows, as expected by regarding the limitations of the structural model used (elastic beam to column connections and no rebar slip), that not all structural stiffness degradation mechanisms are always well represented by the inelastic structural model. In particular, the joints and supports often need special attention and more elaborate models.

For damping model C 1 (t), the good results shown in Fig. 7 have been obtained with α and β computed so that ξ 1 = ξ 2 = 3.3%. For damping model C 2 (t), curves plotted in Fig. 7 have been obtained with the following parameters identified so as to obtain good match between experimental and simulated responses: conditions are satisfied. First, the hysteretic energy which is dissipated by the inelastic response of the structural model has to reach -for the frame considered in this work -150N.m. Then, the seismic demand must be such that the increase in Arias intensity AI [START_REF] Arias | Seismic design for nuclear power plants, chapter A measure of earthquake intensity[END_REF] within the time range [t 1 ; t 1 +10× T ini 1 ] is larger than 0.0025g.s; T ini 1 is the fundamental period of the structure after dead load has been applied and before the earthquake (T ini 1 = 0.36s for the structure considered here).

• Phase 1: from 0 ≤ t ≤ t 1 , ξ 1 = ξ 2 = 1.0%. t 1 is defined such that two
• Phase 2: from t 1 ≤ t ≤ t 2 , ξ 1 = ξ 2 = 4.0%. t 2 is defined as t 2 = t 1 + 10 ×T ini 1 .
• Phase 3: from t 2 ≤ t ≤ T , ξ 1 = ξ 2 = 2.5%, where T is the duration of the seismic signal.

It might happen that one of these two criteria is never satisfied within the earthquake duration. In such a case, only phase 1 is effectively active throughout the analysis.

In the rest of this chapter, we focus only on the new model with changing coefficients, for it is very likely to deliver better prediction from the standard Rayleigh damping models given its more sound physical basis. The damping parameters defined above are expected to be suitable to model the response of the frame structure for every seismic signal that will be used in the next sections, as we usually proceed with damping model C 1 (t). In the next sections, we then investigate likely consequences of using damping model C 2 (t) of C 1 (t) when it comes to structural fragility analysis.

Selection of a set of real ground motions

Following [START_REF] Bommer | The use of real earthquake accelerograms as input to dynamic analysis[END_REF], there are three basic options available for obtaining accelerograms for inelastic time-history analysis: i) to use spectrum-compatible synthetic accelerograms with realistic energy, duration and frequency content; ii) to use synthetic accelerograms generated from seismological source model and accounting for path and site effects; and iii) to use real accelerograms recorded during earthquakes. We chose this latter option and present in the following how we proceed to select ground motion time histories in the PEER ground motion database [START_REF]Users manual for the PEER ground motion database web application[END_REF].

Likely earthquake scenarios in the Cascadia subduction zone

Ground motion time history recorded from the Olympia, Washington 1949 earthquake has been used for the shaking table test at École Polytechnique of Montreal [START_REF] Filiatrault | Seismic performance of ductile and nominally ductile reinforced concrete moment resisting frames. i. experimental study[END_REF]. It is then assumed that the frame structure considered in this investigation is located in the Cascadia subduction zone. The Western margin of the North American plate -from the North of California state up to Vancouver island, British Columbia, Canada -is characterized by the so-called Cascadia subduction zone, where it is subducted by the Juan de Fuca plate beneath the Pacific ocean.

The seismic activity in the Cascadia subduction zone has been investigated for several decades [START_REF] Baker | Source parameters of the 1949 magnitude 7.1 South Puget Sound, Washington, eathquake as determinated from long-period body waves and strong ground motion[END_REF][START_REF] Walter | S. geological survey professional paper 1560 Engineering characterization of earthquake strong motions in the Pacific Northwest[END_REF][START_REF]Damaging capacity of Cascadia subduction earthquakes compared with Chilean subduction[END_REF][START_REF] Katy R Wiest | Source investigation and comparison of the 1939, 1946, 1949 and 1965 earthquakes, cascadia subduction zone, western washington[END_REF]. Three types of earthquakes are produced in this zone:

• Shallow crustal earthquakes are associated to surface faults in the American continental plate with magnitude M w larger than 7.0 and hypocenter depth less than 30 km;

• Thrust interplate or interface earthquakes are due to differential motion in the interface between the Juan de Fuca and the North American plates. They happen offshore with surface hypocenter, generally with depth less than 30 km. The Cascadia subduction zone has the potential to produce a large event of M w = 8.3 ± 0.5;

• Intraplate or intraslab earthquakes occur deep within the Cascadia subduction zone (depth > 40 km) beneath the Puget Sound of Western Washington state. These types of earthquakes have occurred frequently including in 1949 Olympia (M w = 6.8), 1965 Seattle-Tacoma (M w = 6.8), 1999 Satsop (M w = 5.9) and 2001 Nisqually (M w = 6.8);

• No seismicity has been observed for depth larger than 100 km.

Search for a real ground motions dataset

To that purpose, we use the PEER ground motion database with its Web application [START_REF]Users manual for the PEER ground motion database web application[END_REF]. Currently, this database is limited to recorded time series from shallow crustal earthquakes only. A basic criterion used by the Web application to select a representative acceleration time history is that its elastic response spectrum provides a good match to a user target spectrum over a range of periods of interest. We define the target spectrum as the elastic response spectrum corresponding to the feedback accelerogram recorded on the shaking table during the test (see Fig. 2). The Web application allows for assigning different weights to different period ranges so that the matching process is guided by the period ranges with the higher weights. In our case, the response of the structure is governed by its first eigenmode and the fundamental period was experimentally measured to vary within the range 0.36s ≤ T 1 ≤ 0.55s corresponding to the pre-and post-seismic states. We consequently favor the target and actual spectra to match within this period range.

Figure 8: Cross section of the Cascadia subduction zone (adapted from [START_REF] Walter | S. geological survey professional paper 1560 Engineering characterization of earthquake strong motions in the Pacific Northwest[END_REF]).

Table 1: Geologic profile for the Olympic Highway Test Lab strong-motion recording site in Washington [START_REF] Walter | S. geological survey professional paper 1560 Engineering characterization of earthquake strong motions in the Pacific Northwest[END_REF]. v s is the shear wave velocity and ρ the density. Then, several parameters -or acceptance criteria -can be input into the Web application to characterize the likely earthquake scenarios of interest for the study. According to i) what was stated in the previous section 4.1, ii) the cross section of the Cascadia subduction zone depicted in Fig. 8, and iii) the geological profile of the Puget Sound region described in Tab. 1, we define suitable ranges for these parameters as summarized in Tab. 2. Practically, due to the above-mentioned current limitation of the PEER database, we only look for ground motion records corresponding to shallow crustal earthquakes. Table 2: Acceptance criteria for initial ground motions search with the PEER database Web application [START_REF]Users manual for the PEER ground motion database web application[END_REF]. M w is the moment magnitude, R JB is the Joyner-Boore distance, R rup is the closest distance to rupture plane and v s30 is the average shear wave velocity in the top 30 meters of the site. All fault types are considered. Because a first search with these criteria led to a selection of around 20 earthquakes only, we proceeded to a second search with the range of allowed moment magnitudes extended to 6.0 ≤ M w ≤ 9.0. Then, we only retained 48 earthquakes which had the best fitting coefficients with the target response spectrum and finally multiply each of them by a scale factor of 5.0. Note that PEER ground motion database provides the fault normal (FN) and fault parallel (FP) components of the seismic signal and that the 48 records we selected either corresponds to the FN or FP component pertaining to 48 different earthquakes.

Depth [m] Geology Description v s [m/s]

Seismic fragility analysis

In this section, we focus on the vulnerability analysis of the test RC frame structure presented in section 2. Uncertainty is only considered in the seismic loading: the selected time-history ground motions that are likely to occur in the Cascadia subduction zone are used as inputs of inelastic time-history deterministic simulations to compute fragility curves. The damping model -either C 1 (t) or C 2 (t) -is the only variable considered in the RC frame model.

Theoretical background

Structural vulnerability analysis is evaluated here by computing fragility curves which provide the conditional probability for a structure to sustain a given damage level for a given earthquake intensity. Following [START_REF] Masanobu Shinozuka | Statistical analysis of fragility curves[END_REF], it is assumed that the fragility curve can be expressed in the form of a two-parameter lognormal distribution function. The estimation of these two parameters is then performed with the maximum likelihood method.

Let consider a sample of n independent observations x j that can be classified in two classes as "success" (x j = 1) and "failure" (x j = 0). Let also consider that each realization x j has a relative frequency of success equal to F (IM j , θ), that is a function that depends on the intensity measure IM j associated to each realization x j and of a set of parameters θ. Then, the probability of observing a set of realizations x = (x 1 , ..., x n ) composed of p successes x 1 , ..., x p and n -p failures x p+1 , ..., x n , whatever the order, is ( [START_REF] Poincaré | Calcul des Probabilités[END_REF], p.77)1 :

f (x|θ) = n! p! (n -p)! (F (IM j , θ)) p • (1 -F (IM j , θ)) n-p (10)
The problem one has to solve can be expressed as: Given the observed data x, find the set of parameters θ that is most likely to have produced these observed data. To solve this inverse problem, we define the likelihood function as a function of θ given x:

L(θ|x) = f (x|θ). ( 11
)
The principle of maximum likelihood estimation states that: given the data x actually observed, the set of parameters θ looked for is the one that makes x the most likely data to be observed. θ can thus be identified by maximizing the likelihood function L. For computational convenience, the log-likelihood lnL(θ|x) is introduced and maximized, which provides the same estimators because ln x is a monotonic function. The problem one has to solve thus reads:

θ e = argmax θ lnL(θ|x) (12)
Under the lognormal assumption, the fragility curve for a particular damage index DI i -defining what is "success" and "failure" -is defined as:

F (IM j , θ e ) = φ 1 ζ e ln IM j c e = P [DI ≥ DI i |IM j ] (13)
where θ e = {ζ e , c e } is the set of estimated parameters and φ(•) is the standardized normal distribution function.

Constructing fragility curves in such a framework raises issues concerning their statistical significance. In [START_REF] Masanobu Shinozuka | Statistical analysis of fragility curves[END_REF], Shinozuka et al. provide tools to test the goodness of fit between the inferred fragility curve and the realization of the random variable X j following Bernoulli distribution: X j = 1 when the damage index is reached and X j = 0 otherwise. They also present a Monte Carlo technique they use to demonstrate the extent of the statistical variations in the estimators θ. Another very important contribution for the statistical analysis of fragility curves is the work of Sáez et al. [START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF]. They show how to compute the amount of Fisher information about the set of parameters θ -the terms of the Fisher information matrix being by definition F ij (θ) = cov ∂lnL ∂θ i ; ∂lnL ∂θ j

, where cov(•) denotes the covariance -which is provided by the ground motions used to construct fragility curves, and then how to compute a lower bound for the standard deviation of the elements of θ e . This method thus provides a way to somehow measure the ability of the data to estimate θ.

Intensity measures and structural damage indices

To characterize the seismic ground motion time-histories used for structural vulnerability analysis, they are assigned intensity measures (IM). 18 of them are reviewed in [START_REF] Léger | Incremental dynamic analysis of nonlinear structures: Selection of input ground motions[END_REF] in the context of the issue of selecting earthquakes for incremental dynamic analysis of inelastic steel frame structures. 44 IM are reviewed or proposed for masonry structures in [START_REF] Bommer | The influence of strong-motion duration on the seismic response of masonry structures[END_REF]. Among the most common IM, one finds: the peak ground acceleration (P GA); the acceleration at the fundamental period S A (T 1 , 5%); the Arias intensity (AI) which is effectively a measure of the total energy in the ground motion and computed as

AI(t) = π 2g t 0 ü2 g (τ )dτ ; (14)
the significant duration D 5-95 , which the time needed to build up between 5% and 95% of the total Arias intensity. Other measures are: the root mean square of acceleration (RMSA) computed as [START_REF] Léger | Incremental dynamic analysis of nonlinear structures: Selection of input ground motions[END_REF]:

RMSA = 1 τ d 2g π AI, (15)
where τ d is an effective duration of the record, taken here as D 5-95 ; spectrum intensities such as [START_REF] Bommer | The influence of strong-motion duration on the seismic response of masonry structures[END_REF]:

SI V = T b Ta S V (T, ξ) dT, (16)
where S V is the spectral velocity at period T and with damping ξ. Here, we take T A = 0.36s and T B = 0.55s which are the fundamental periods in pre-and postearthquake conditions. For fragility surfaces, Seyedi et al. [START_REF] Seyedi | Development of seismic fragility surfaces for reinforced concrete buildings by means of nonlinear time-history analysis[END_REF] chose the spectral displacement S D (T 1 ) and S D (T 2 ) where T 1 and T 2 are the two main eigenperiods in the direction along which the ground motions are applied.

Then, to characterize the structural response, damage indices (DI) are used. One can distinguish between three categories: i) DI computed from energy quantities such as the Park-Ang-Wen damage accumulation model [START_REF] Park | Seismic damage analysis of reinforced concrete buildings[END_REF] or the normalized hysteretic energy used in [START_REF] Castiglioni | Failure criteria and cumulative damage models for steel components under cyclic loading[END_REF] or [START_REF] Léger | Incremental dynamic analysis of nonlinear structures: Selection of input ground motions[END_REF]; ii) Other DI based on quantities directly related to the structural inelastic response such as the ductility demand [START_REF] Léger | Incremental dynamic analysis of nonlinear structures: Selection of input ground motions[END_REF], the cumulative ductility index as defined in [START_REF] Castiglioni | Failure criteria and cumulative damage models for steel components under cyclic loading[END_REF] or the maximum strength degradation ratio [START_REF] Bommer | The influence of strong-motion duration on the seismic response of masonry structures[END_REF]; iii) Other quantities not necessarily pertaining to -but affected by -the inelastic structural behavior such as the maximum relative roof displacement (MRD) or the maximum inter-story drift ratio (MISDR). To these two latter indices, because particular attention is paid in this work on the energy dissipated by both the damping and the inelastic structural models, we also compute the amount of hysteretic energy (E H ) dissipated by the inelastic mechanisms explicitly accounted for in the structural model and the amount of damping energy (E D ), as well as their respective ratio E H,D R = E H,D /(E H + E D ).

Fragility curves

We first investigate in Table 3 the correlation between the intensity measures and damage indices considered in this work, when a linear model is used to predict a DI from an IM. Correlation is the lowest for intensity measure D 5-95 because there is no explicit influence of the earthquake duration in the various DI considered. The correlation for the energy ratios is poor for all the IM considered, which is not the case for the dissipated energy quantities where, in particular, correlation is good with the energy contained in the seismic signal (AI). The very good correlated maximum roof displacement and spectral displacement at fundamental period T 1 with 5% viscous damping -as depicted in Fig. 9 -will be used in the following to construct fragility curves. In Fig. 9, one can see an obvious tendency to obtain larger maximum roof displacement with damping model C 2 (t). Constructing fragility curves for a given MRD level along with a proper statistical analysis provides another way to infer some likely consequences of using a damping model instead of another. To that purpose, Fig. 10 shows the fragility curves (solid lines) pertaining to DI MRD ≥ 40mm obtained with both damping models C 1 (t) and C 2 (t). Dashed lines approximately represent the smallest area the fragility curves would describe when its parameters θ = {ζ, c} varies around θ e = {ζ e , c e } plus or minus their standard deviation. Dashed lines are constructed as in [START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF], relying on the inverse of the Fisher information matrix provided by the selected ground motion about θ.

From Fig. 10, one can infer that the additional damping model entails very

Conclusions and perspectives

The inelastic structural response in seismic loading results from the combination of an inelastic structural model with an added damping model. The damping model has to be consistent with the inelastic structural model, which implies that it has, and only has, to model the energy dissipation sources not otherwise explicitly accounted for in the inelastic structural model, nothing more or less. Rayleigh damping models are the most commonly used for earthquake engineering applications, although it is well established that controlling the amount of energy dissipation these models introduce throughout inelastic time history simulations is difficult to achieve. On another hand, inelastic time history analyses are widely used for structural fragility assessment in seismic loading. Uncertainties arising from the seismic signal likely to excite a building in a given geographical region as well as from the lack of knowledge on the structural parameters which characterize the inelastic structural model have both been considered in structural fragility analyses. On the opposite, the damping model is scarcely considered as a source of uncertainty; moreover, when it is the case, the damping ratio is assigned a probabilistic distribution but whether commonly used added damping models are suitable or not is an issue that is not regarded. This is the issue addressed in the work presented above.

In this chapter, we propose a new family of Rayleigh damping models that relies both on the capacity of the inelastic structural model to absorb energy and on the seismic demand and compare it to a commonly used Rayleigh model in the context of the fragility analysis of a RC moment-resisting frame. Its inelastic response is modeled by fiber elements using a constitutive law recently developed by the authors [START_REF] Jehel | Towards robust viscoelastic-plastic-damage material model with different hardenings / softenings capable of representing salient phenomena in seismic loading applications[END_REF]. The structure is supposed to be built in the Cascadia subduction zone, a seismically active zone in the South-West of Canada and the North-West of the USA. From this comparative analysis, it can be inferred that the added damping model entails very significant uncertainty in structural fragility analysis of inelastic RC frame structures. and C 2 (t). Dashed lines are plotted according to the method presented in [START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF] and show that, for the DI level considered here, the fragility curves are constructed from a number of ground motions which is sufficient for drawing pertinent conclusions from this comparative study. Circles represent the realizations of the Bernoulli random variable X i = 1 when RMD i ≥ 40mm and X i = 0 otherwise, i ∈ [1; 48].

As a further development, we seek the proposed new family of Rayleigh damping models to be confronted to other experimental evidence, criticized and improved. Albeit constructed on a stronger physical basis than commonly used damping models, there is yet no clear guarantee that this model remains consistent with the inelastic structural model it is coupled to, throughout inelastic time history analysis. Finally, such a proposition for a damping model should not eclipse the need for further improvement of the inelastic structural model, which would, for the case treated in this work, at least involve explicitly accounting for the inelastic response of the beam to column connections. ing table tests used in this work. The first author benefited from partial funding from Électricité de France (EDF) within the research project "MARS" ("Modèles Avancés pour le Risque Sismique").

Figure 1 :

 1 Figure 1: RC frame structure tested on the shaking table at École Polytechnique in Montreal. Dimensions are in [mm].

1 and T ini 1 coincideFigure 4 :

 14 Figure 4: Internal variables for compression and their phenomenological interpretation. [top, left] Elastic response until yield stress σ y is reached. [top, centre] In the plastic domain, plastic deformation (ǭ p ) and strain hardening ( ξp ) develop. [top, right] Once the limit stress σ f is reached, damage mechanisms are activated too, leading to a degradation of stiffness ( D) and a change in the strain hardening evolution ( ξd ). [bottom, left] Once the ultimate stress σ u is reached, deformation localizes ( ūp ) and strain-softening is observed ( ξp ). [bottom, center] Local hysteresis loops are represented with kinematic hardening in the plastic domain ( λp ). Note that each parameter of this model has a clear interpretation regarding the constitutive law to be identified.

Figure 5 :

 5 Figure 5: Finite element mesh (dimensions in [mm]) and material constitutive laws for the inelastic structural modeling. [bottom left] Confined and unconfined concrete behavior laws. [bottom right] Steel constitutive law.

Figure 6 :

 6 Figure 6: Three common phases in the seismic response and the concept of "key window". [top left] Structural relative displacement time-history. [top right] Total relative seismic energy time-history in the structure. [bottom] Arias intensity of the seismic signal: AI(t) = π 2g t 0 ü2 g (τ )dτ .

Figure 7 :

 7 Figure 7: Experimental and simulated top-displacement; simulated hysteretic plus elastic energies (E H ), simulated additional damping energy (E D ), and both experimental and simulated total internal work. [top] With common added damping model C 1 (t); [bottom] With the proposed new family of Rayleigh damping models C 2 (t). The structural responses shown here for damping models C 1 (t) and C 2 (t) looks very similar because both models have been calibrated to experimental data; however, model C 2 (t) has more capability for representing transient evolution of added damping.

  [ 0 , 150 ] [ 30 , 200 ] [ 30 , 200 ] R rup [km] [ 0 , 150 ] [ 30 , 100 ] [ 30 , 100 ] v s30 [m/s] [ 0 , 200 ] [ 0 , 200 ] [ 0 , 200 ]

Figure 9 :

 9 Figure9: Linear model between S A (T 1 , 5%) and MRD for both damping models C 1 (t) and C 2 (t).

1 S

 1 A (T 1 , 5%) [g]

Figure 10 :

 10 Figure10: Structural fragility curves (solid lines) for both damping models C 1 (t) and C 2 (t). Dashed lines are plotted according to the method presented in[START_REF] Patricio | Effect of the inelastic dynamic soil-structure interaction on the seismic vulnerability assessment[END_REF] and show that, for the DI level considered here, the fragility curves are constructed from a number of ground motions which is sufficient for drawing pertinent conclusions from this comparative study. Circles represent the realizations of the Bernoulli random variable X i = 1 when RMD i ≥ 40mm and X i = 0 otherwise, i ∈ [1; 48].

Figure 11 :

 11 Figure 11: [left] Analysis where only phase 1 is activated for damping model C 2 (t). [right] Analysis where three phases are accounted for in damping model C 2 (t). [top] Ground motion time-history. [bottom] Top relative displacement time-history; Model C 1 (t) with ξ 1 = ξ 2 = 3.3% all along the simulation; [bottom left] Model C 2 (t) with ξ 1 = ξ 2 = 1% all along the simulation; [bottom right] Model C 2 (t) with ξ 1 = ξ 2 = 1% for 0 ≤ t ≤ 15.4s, ξ 1 = ξ 2 = 4% for 15.4s ≤ t ≤ 19s and ξ 1 = ξ 2 = 2.5% for 19s ≤ t ≤ 60s.

Table 3 :

 3 Square of the coefficient of correlation R 2 (in %) between actual DI and DI obtained from a linear model between IM and DI, for added damping models C 1 (t) / C 2 (t). Fig.11is shown to provide better insight into the discrepancies one can expect to observe in the structural response when either damping model C 1 (t) or C 2 (t) is used. Both seismic signals considered have the same intensity measure S A (T 1 , 5%) = 0.51g. For the analysis in concern on the left part of Fig.11, either the maximum capacity of the structure to store hysteretic energy E H = 150N.m or the minimum seismic demand 0.0025/(t 2 -t 1 ) is not reached, so that, for damping model C 2 (t), only the phase 1 with ξ 1 = ξ 2 = 1% is effectively active throughout the analysis; this makes both damping models analogous but with critical damping ratios in sharp contrast -3.3% against 1% -leading to very different structural responses. For the analysis in concern on the right part of Fig.11, phase 2 begins at t 1 = 15.4s until t 2 = 19s. Both models predict almost identical MRD but oscillation amplitude is more rapidly attenuated for model C 1 (t).

		MRD MISDR E D/H R	E D	E H
	P GA	73 / 71 71 / 69 3.5 / 0.6 51 / 52 36 / 37
	AI	68 / 68 68 / 68 9.1 / 5.9 86 / 83 61 / 62
	D 5-95	8 / 7	8 / 7	1.3 / 1.0	1 / 1	2 / 2
	RMSA	75 / 73 74 / 73 8.2 / 4.5 57 / 55 43 / 45
	S A,V,D (T 1 , 5%) 85 / 82 84 / 82 5.7 / 2.6 68 / 68 53 / 55
	SI A,V,D (5%)	77 / 74 77 / 75 8.5 / 4.6 57 / 57 48 / 51
	significant uncertainty in structural fragility analysis of inelastic RC frame struc-
	tures.

Henri Poincaré (1854[START_REF] Poincaré | Calcul des Probabilités[END_REF]) is a French mathematician, physician and philosopher. This year is the hundredth anniversary of his death.

Acknowledgements

The authors thank Pr. André Filiatrault for providing the data from the shak-


Christophe Kervazo
email: christophe.kervazo@cea.fr

Jérôme Bobin

Cécile Chenot

Heuristics for Efficient Sparse Blind Source Separation

Sparse Blind Source Separation (sparse BSS) is a key method to analyze multichannel data in fields ranging from medical imaging to astrophysics. However, since it relies on seeking the solution of a non-convex penalized matrix factorization problem, its performances largely depend on the optimization strategy. In this context, Proximal Alternating Linearized Minimization (PALM) has become a standard algorithm which, despite its theoretical grounding, generally provides poor practical separation results. In this work, we first investigate the origins of these limitations, which are shown to take their roots in the sensitivity to both the initialization and the regularization parameter choice. As an alternative, we propose a novel strategy that combines a heuristic approach with PALM. We show its relevance on realistic astrophysical data.

I. INTRODUCTION

A. Blind Source Separation problem

In the BSS [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF] framework, the data are composed of m observations, each of which has t samples. These observations are supposed to be some linear combinations of n sources. The objective of BSS is to retrieve the sources as well as the mixing coefficients. In matrix form, the goal is therefore to find two matrices S (of size n × t) and A (of size m×n), called respectively the source and the mixing matrices, such that: X = AS + N, where X (of size m × t) is the observation matrix that is corrupted with some unkwown noise N. Since it requires tackling an ill-posed unsupervised matrix factorization problem, further assumptions are needed, including the statistical independance of the sources (ICA - [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]), the non-negativity of A and S [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF]. In this work, we will focus on the sparsity of the sources [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF]- [START_REF] Li | Underdetermined blind source separation based on sparse representation[END_REF]. In this framework, sparse BSS will aim at finding a (local) minimum of:

min A,S 1 2 X -AS 2 F + ι {∀i; A i 2 2 ≤1} (A) + R S (SΦ T S ) 1 (1)
The first data fidelity term promotes a faithful reconstruction of the data. The use of the Frobenius norm . F stems from the assumption of a Gaussian noise N. The second term involving the characteristic function ι corresponds to the oblique constraint ensuring that the columns A i of A are all in the 2 ball. This avoids degenerated A and S matrices where A → ∞ and S → 0. The last term involving the Hadamard product enforces a 1 sparsity constraint in a transformed domain Φ S . In the following, Φ S will be taken equal to either the identity (in which case the sparsity is enforced in the direct domain) or the starlet transform [START_REF] Starck | Sparse Image and Signal Processing: Wavelets, Curvelets, Morphological Diversity[END_REF]. R S controls the trade-off between the data fidelity and the sparsity terms. It can be decomposed into R S = Λ S W where Λ S is a diagonal matrix of the regularization parameters λ 1 , λ 2 , ..., λ n and W is a matrix used to introduce individual penalization coefficients in the context of reweighted 1 [START_REF] Candes | Enhancing sparsity by reweighted 1 minimization[END_REF] (when no reweighting is used, W is equal to the identity matrix).

B. Sparse BSS in practice

Since sparse BSS requires solving a penalized matrix factorization problem, it is important to highlight that the separation quality strongly depends on the optimization strategy. For that purpose, different algorithmic frameworks have been used so far: projected Alternate Least-Square (ALS - [START_REF] Gillis | Accelerated multiplicative updates and hierarchical ALS algorithms for nonnegative matrix factorization[END_REF]), Block-Coordinate Descent (BCD - [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]) and PALM [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. However, any practitioner can draw the same conclusion: the solution of sparse BSS methods is highly sensitive to the initial point and the values of the regularization parameters, which are generally tricky to tune without any first guess of the solution.

Initialization: As problem (1) is not convex but multiconvex, the algorithm performing its minimization can be trapped in spurious critical points, depending on the initial matrices A and S. Different optimization strategies can be affected differently.

Regularization parameters: The practical choice of the regularization parameters Λ S is of paramount significance:

• The parameters are directly impacting the shape of the solution through the trade-off between the sparsity level and fidelity to the data. • The problem being non-convex, a change in the parameters can also bring the optimization algorithm to stabilize in the neighborhood of a different critical point, potentially leading to a very different solution.

C. Contributions

While PALM is theoretically well rooted and yields rather fast minimization schemes (in contrast to BCD), it generally provides poor separation results. We first investigate the origins of this poor practical efficiency, and especially the sensitivity to both initialization and regularization parameter settings. We further show how PALM-based implementations can benefit from the information provided by heuristic approaches which are in pratice more robust. The robustness of the proposed combined strategy is demonstrated on realistic astrophysical data.

II. OPTIMIZATION STRATEGY: PALM AND GMCA

As emphasized before, the optimization strategy is decisive to avoid stationary points. In the following, we will discuss PALM [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and projected ALS such as Generalized Morphological Component Analysis (GMCA - [START_REF] Bobin | Blind source separation: The sparsity revolution[END_REF]). BCD will not be studied since it requires an exact minimization at each iteration, leading in general to a high computational cost.

A. PALM algorithm PALM has been introduced in [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF] and has benefited from several extensions [START_REF] Chouzenoux | Variable metric forwardbackward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF], [START_REF]A block coordinate variable metric forward-backward algorithm[END_REF]. Its attractiveness partially comes from the fact that it is proved to converge to a local minimum of (1) under mild conditions [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. PALM is an iterative algorithm. At each iteration, it alternates between a proximal gradient step on A and S. In our case, it can be shown that the proximal operator corresponding to the 1 sparsity term in (1) is the soft-thresholding operator (with for each source S λ (.) = sign(.) [|.| -λ] + ) applied in the transform domain. Furthermore, the operator corresponding to the oblique constraint is the projection of each column of A on the 2 unit ball, which we shall denote as Π . 2 (.). Each iteration (k) can then be read as:

1 -Update of S using a fixed A:

Ŝ = S (k-1) - γ A (k-1)T A (k-1) 2 A (k-1)T (A (k-1) S (k-1) -X) (2) S (k) = S γR S A (k-1)T A (k-1) 2 ŜΦ T S Φ S (3)
2 -Update of A using a fixed S:

 = A (k-1) - δ S (k) S (k)T 2 (A (k-1) S (k) -X)S (k)T (4) A (k) = Π . 2  (5)

B. GMCA algorithm

GMCA algorithm is based on projected ALS. At each iteration (k), the gradient step of ( 2) and ( 4) is replaced by a multiplication by a pseudo-inverse: 1 -S is updated assuming a fixed A.

Ŝ = A (k-1) † X (6)
Where 1) .

A (k-1) † is the pseudo-inverse of A (k-
S (k) = S R S ( ŜΦ T S )Φ S (7)
2 -A is updated assuming a fixed S:

 = XS (k) † (8)
A (k) = Π . 2 ( Â) (9)
Compared to PALM, GMCA is only a proxy for the minimization of (1) and cannot be proved to converge (while stabilizing in most of practical cases after some iterations). However, thanks to heuristics, it benefits from an automatic thresholding strategy (see Sec.IV-A1) which has been empirically shown to improve its robustness with respect to the initialization.

III. PALM-BASED SPARSE BSS IN PRACTICE

The objective of this section is to empirically shed light on the lack of robustness of PALM to motivate the need of finding a relevant strategy to enable its use on practical problems. The first subsection presents our experimental protocol, while the second and third ones focus on the sensibility of PALM with respect to the regularization parameters and the initialization.

A. Experimental protocol

To bring out PALM sensibility and the mechanisms at stake, our experiments focus in this subsection on simple data coming from n = 2 equilibrated (with equal 2 norm) sources. The sources are assumed to be exactly sparse in the direct domain, a proportion p = 0.1 of the t = 500 samples being nonzeros and drawn according to a standard normal distribution. The mixing is performed through a A matrix drawn randomly following a standard normal distribution and modified to have unit columns. Its condition number is C d = 10. There are m = n observations. To complete the creation of the X data, a Gaussian additive noise is added to the mixing, such that the Signal to Noise Ratio is SNR = 60 dB. Once the algorithm launched, it stops based on a convergence criterion computed through the angular distance of the columns of A between two successive iterations. The separation quality is then mesured using a global criterion on A [START_REF] Bobin | Sparsity and adaptivity for the blind separation of partially correlated sources[END_REF]: C A = -10 log 10 (median(|PA † A * | -I d )), with A * the true mixing matrix and PA † the pseudo-inverse of the solution found by the algorithm corrected through P for the scale and permutation indeterminacies. The higher C A , the better the separation.

To study PALM robustness, we perform an exhaustive search for the parameters, testing many different values. The feasability of this approach is due to both the fact that we work on simulated data, for which we know the true A and S, and the low number of sources, asking for only two parameters. Such an approach would be intractable on practical cases.

B. Sensibility to the thresholding strategy

Following our experimental protocol, we obtain Figure 1, on which C A is displayed as a function of the threshold of each of the 2 sources ((λ 1 , λ 2 ) are part of R S and therefore impact PALM through Eq. 3). On this figure, the lack of robustness and sensibility to the parameter choice is underlined by the following remarks:

C. Sensibility to the initialization

As described in section I-B, the initialization impacts the results due to the non-convex nature of problem (1) because the algorithm can stay trapped in spurious critical points. The standard deviation of C A over different initializations of A is plotted in Fig. 2 to quantify this impact. It can reach up to 5 dB for some regularization parameters choices. The best parameter setting corresponds to a quite high standard deviation and can shift depending on the initialization. It makes the practical choice of an initialization a critical point which directly impact the quality of the local minimum found by the algorithm. It further highlights the sensitivity to the values of the regularization parameters. A key question is thus how to tune PALM in practice and if it can be done in an automatic way without prior knowledge about the solution.

IV. COMPLEXITY OF INTRODUCING HEURISTICS IN PALM

A. Heuristic motivation and description

Building on the automatic thresholding strategy of GMCA, the goal of this part is to try to derive one for PALM.

1) Automatic parameter choice in GMCA: In GMCA, the thresholding is directly performed on a minimizer of an approximation of the data fidelity term. Said differently, it is performed in an approximation of the source space. This has lead to a simple interpretation of the parameters, allowing for an heuristic thresholding strategy that yields good practical results and robustness. The true sources being sparse, the thresholding should remove a dense signal, partially coming from the back-projection in the source space of the data noise N. In practice, this dense signal removal is performed computing the Median Absolute Deviation (MAD) for each source and setting the corresponding threshold to a multiple k of this value:

(λ1, λ2, ..., λn) T = k × MAD( Ŝ) (10)
In this equation and in the remaining of this section, Φ T S was supposed to be the identity matrix without loss of generality. k is a positive number and the MAD operator is computed row-byrow with MAD(z) = median(|zmedian(z)|) for z ∈ R t . If we suppose that GMCA has found the true A matrix denoted by A * , it can be shown that:

k × MAD( Ŝ) k × MAD A * † N (11)
It has however to be emphasized that this data noise removal interpretation of the thresholding is limited to the case when the estimation of the mixing matrix is close to A * . In the opposite case, the imperfect demixing results in an increase of the thresholds due to an extra noise coming from interferences between the sources.

2) Heuristic in PALM: Let us assume that the thresholds in PALM are computed the same way as in GMCA through eq. ( 10), and that the algorithm has converged to both the true matrices A * and S * . Then:

(λ1, λ2, ..., λn) T = k × MAD S * - γ A * T A * 2 A * T (A * S * -X) k × γ A * T A * 2 MAD(A * T N) (12)
where the last line is obtained thanks to the assumption of S * sparsity. Therefore, using the MAD enables a thresholding of a projection of the noise N, which yields a similar interpretation as in GMCA. This observation must however be tempered:

B. Experimental protocol

The mixing and the metric on the separation quality are the same as in Sec. III-A. The automatic parameter setting described in the previous subsection is used inside of a PALM algorithm for different k values. k is the same for both sources, since the MAD is supposed to be adaptive enough to the noise.

C. Empirical performances of the heuristic

The results are displayed in Fig 3 . The low values of C A , reaching at most 6dB, are to be compared with the 39 dB obtained during the exhaustive parameter search in III-A. While efficient in GMCA, the MAD heuristic leads to a bad separation quality due to the different space in which the thresholding is performed (cf. IV-A2), yielding a different role to the thresholds. To better understand this role, let us assume that the algorithm is initialized with the true mixing matrix A * only (and not with S * , contrary to what is described in IV-A2). In PALM, the use of the transpose A * T in eq. ( 2) leads towards a re-mixing of the sources at each iteration whose goal is to gradually decrease the data fidelity term. This re-mixing however increases the sparsity term and should be impeded by the thresholds, giving to them a demixing role. On the opposite, when using GMCA, such a re-mixing between the sources does not exist due to the use of the pseudoinverse A * † instead of A * T , which alleviate somehow the issue of finding relevant thresholds. Therefore, using the MAD heuristic in PALM does not results in good separations in practice, because the interpretation in terms of noise removal is less relevant than in GMCA.

V. COMBINING GMCA AND PALM: A HYBRID STRATEGY

In this part, we propose to combine the best of GMCA and PALM in a two step approach. The algorithm consists in a warm-up stage, in which GMCA is performed, followed by a refinement stage during which PALM is performed retaining as much information as possible coming from the warm-up stage.

A. Motivation and full description of the algorithm

Our approach is motivated by several remarks:

• PALM theoretical background: while GMCA is only a proxy, the 2-step algorithm will attempt to solve exactly eq. ( 1) as PALM does. It further benefit from the high potential accuracy of PALM. Moreover, PALM is proved to converge under mild assumptions [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], so is the 2-step algorithm.

• GMCA robustness with regards to initialization: we propose to re-use GMCA outputs A GMCA and S GMCA as an initialization of PALM. Due to the almost constant results of GMCA when confronted to different initializations, the quality of the output of PALM in the refinement stage will have a much decreased variance. • Benefit from GMCA solution: the results of GMCA can be re-used by PALM in the refinement stage as a first guess about the shape of the solution. It is done through two mechanisms: -Since both A GMCA and S GMCA are close to A * and S * , they can be used to derive the thresholds using the MAD. It has however to be emphasized that after a given number of iterations the thresholds become fixed, which is a necessary condition to ensure convergence. -Furthermore, S GMCA should already give a good approximation of the most prominent peaks. This information can be exploited in the refinement stage through the introduction of reweighted L1 [START_REF] Candes | Enhancing sparsity by reweighted 1 minimization[END_REF] using the reweighting matrix W in problem [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]. In brief, the thresholds are lowered for the biggest samples of the estimated sources, reducing the bias introduced by soft-thresholding. This is particularly appealing in PALM because the biases in these coefficients behave as errors that are transmitted to the other sources through the remixing introduced by A T . These remarks lead to the following algorithm using a smart initialization and thresholding strategy:

Input : X (data matrix)

• Random initialization A 0 and S 0 • Warm-up stage: A GMCA , S GMCA = GMCA(X,A 0 ,S 0 ) • Refinement stage:

A PALM , S PALM = PALM(X,A GMCA ,S GMCA ) The initialization, thresholding strategy and reweighting information come from the warm-up stage.

B. Empirical results

The experimental protocol is the same as described in Sec. IV-B except that the 2-step algorithm is used instead of PALM. The results are plotted in Fig. 3. With values of C A higher than 33 dB, the demixing is close to the best ones obtained with the exhaustive search in Sec. III. Compared to PALM only, the variance of results over different initializations is also much decreased, becoming almost zero, which shows the robustness of the algorithm with regards to the initialization.

VI. EXPERIMENT ON REALISTIC DATA

The goal of this part is to apply our algorithm to realistic data to show its efficiency.

A. Data description and experimental protocol

The n = 2 sources come from simulations obtained from real data of Cassiopeia A supernova remnant. These data originate from the Chandra X-ray observatory. The sources in these wavelength values correspond to the thermal emission and the iron. As displayed in Fig. 4, they each consists in a 2D image of resolution t = 128 × 128 pixels, supposed to be approximately sparse in the starlet domain. The mixing matrix is simulated as described in III-A. To increase the realism of the data and further test the algorithm, we tried three relatively low SNR values: 10, 15 and 20 dB. The sparsity is enforced in the starlet domain. k is set to 3, which corresponds to a standard hypothesis in terms of Gaussian noise removal.

B. Empirical results

The mean of C A values is displayed in Fig. I. To visually assess the separation quality, the estimation for a SNR of 15 dB is shown in Fig. 4. The 2-step approach always achieve better results than the two classical BSS algorithms with which we performed the comparison, namely Relative Newton Algorithm (RNA) and Efficient FastICA (EFICA). It also outerpeforms both GMCA and PALM, being always better by at least 2 dB than the best of them. The bigger differences with GMCA than in Fig. 3 are probably explained by the fact that the approximation made by GMCA is less well verified due to both a higher noise level and less sparse sources. In addition to the results displayed in Fig. (I), the standard deviation of C A over different initializations is almost 0, which shows the robustness of the algorithm.

CONCLUSION

In this work, we show that the ability of recent and theoretically grounded optimization strategies like PALM to provide solutions to sparse BSS problems is highly sensitive to both the initialization and the values of the regularization parameters. To further design an effective as well as robust algorithm, we introduce a 2-step strategy combining PALM with robust heuristic methods such as GMCA. Beyond improving the robustness of PALM-based implementations with respect to initialization, the regularization parameters can be automatically set in the proposed approach. Numerical experiments on both simulated and realistic data demonstrate a high separation quality and good robustness on mixings with low SNR.

Fig. 1 .

 1 Fig. 1. Mean of C A as a 2D function of the thresholds for 5 different initializations of PALM algorithm.

Fig. 2 .

 2 Fig. 2. Standard deviation of C A as a 2D function of the thresholds for 5 different initializations of PALM algorithm. and λ 2 = 3.88 × 10 -3 . Merely setting λ 1 = λ 2 = 3.81 × 10 -3 yields a loss of more than 10 dB. Second, in this simple toy example, the sources are equilibrated. It would not necessarily be the case with real-life sources, which could lead to a shift of the diagonal; • Inside of the diagonal itself, a 2×10 -3 change of λ 1 = λ 2 can lead to a loss of about 7 dB.

Fig. 3 .

 3 Fig. 3. Mean and quartiles of C A for 10 different initializations of 3 different algorithms : PALM, GMCA and 2 steps.

Fig. 4 .

 4 Fig. 4. Up: true sources. Down: sources estimated by the two step algorithm.

TABLE I AVERAGE

 I C A (10 DIFFERENT INITIALIZATIONS) FOR 3 SNR VALUES AND 5 ALGORITHMS. THE BEST VALUES ARE IN BOLD.

		10 dB	15dB	20 dB
	2 step	11.57	16.92	22.09
	PALM	9.38	10.94	11.01
	GMCA	8.87	12.15	15.87
	EFICA	-6.92	5.11	9.41
	RNA	-6.68	-5.49	6.27

• High quality separation happens for values close to the diagonal. A slight deviation from it can cause a 30 dB drop for a 3 × 10 -4 modification of one parameter and lead to a very bad separation. One could argue that this problem could be circumvented by imposingλ 1 = λ 2 .First, it has to be emphasized that the maximum is not on the diagonal but slightly outside, for λ 1 = 3.81×10 -3

• It only holds when and if PALM has converged towards good A and S; • The projection is not performed in the same space as GMCA, for which the noise is back-projected in the source space. Both spaces however merge when A is orthogonal.

ACKNOWLEDGMENT This work is supported by the European Community through the grant LENA (ERC StG -contract no. 678282), and by the MOD University Defence Research Collaboration (UDRC) for signal processing in the networked battlespace, number P/K014277/1.