Enora Denimal
email: enora.denimal@inria.fr

Jean-Jacques Sinou
email: jean-jacques.sinou@ec-lyon.fr

Samuel Nacivet
email: samuel.nacivet@stellantis.com

Prediction of Squeal Instabilities of a Finite Element Model Automotive Brake With Uncertain Structural and Environmental Parameters With a Hybrid Surrogate Model

Keywords: hybrid surrogate model, stability, friction-induced vibrations, uncertainties

This study focuses on the prediction of the stability behavior of an industrial automotive brake system under structural and environmental uncertainties. Uncertainties are modeled with a random distribution or an interval and are propagated with a hybrid surrogate model associating polynomial chaos and kriging. The objective is to create a surrogate model of each eigenvalue computed with the complex eigenvalue analysis (CEA). As the modes can be tracked only when unstable, the effective size of the training sets can become extremely small. Despite this limitation, it is shown the hybrid meta-model is still able to predict the stability of the brake system. Moreover, the hybrid meta-model gives a direct access to the mean and variance of the eigenvalues with respect to the design parameters without any additional Monte Carlo simulations (MCS). By considering different probability density function for the friction coefficient, it is shown it has a high influence on the stability and the latter should be accurately estimated.

Introduction

Friction-induced vibrations are a complex and rich topic of research since several decades for both industrials and academics [START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part I: Mechanics of Contact and Friction[END_REF][START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part II: Dynamics and Modeling[END_REF]. It is well-known that the phenomenon is fugitive and depends on numerous parameters, such as the lining materials [START_REF] Eriksson | Surface Characterisation of Brake Pads After Running Under Silent and Squealing Conditions[END_REF], the wear [START_REF] Bergman | Influence of Disc Topography on Generation of Brake Squeal[END_REF], the environment [START_REF] Eriksson | A Study of the Influence of Humidity on the Friction and Squeal Generation of Automotive Brake Pads[END_REF], the geometry of the contact surface [START_REF] Bergman | The Effect of Reduced Contact Area on the Occurrence of Disc Brake Squeals for An Automotive Brake Pad[END_REF], etc., making the study and the modeling complex. Thus, the consideration of the variability or the uncertainty associated to some parameters emerged a few years ago. Indeed, in 2009, Culla and Massi [START_REF] Culla | Uncertainty Model for Contact Instability Prediction[END_REF] considered some parameters as uncertain and performed Monte Carlo Samplings (MCS) on a numerical model to improve the squeal prediction. In 2010, Butlin and Woodhouse [START_REF] Butlin | Friction-Induced Vibration: Quantifying Sensitivity and Uncertainty[END_REF] investigated the variability of squeal related parameters and laid the foundations of the consideration of uncertainty for improving squeal predictions. In 2014, Tison et al. [START_REF] Tison | Improvement in the Predictivity of Squeal Simulations: Uncertainty and Robustness[END_REF] demonstrated that considering uncertainty improves the squeal prediction and confirmed their results with experiments. Since, numerous studies were dedicated to the consideration and the modeling of uncertainties in brake systems from a numerical and/or experimental point of view [START_REF] Lü | Brake Squeal Reduction of Vehicle Disc Brake System With Interval Parameters by Uncertain Optimization[END_REF][START_REF] Nobari | Uncertainty Quantification of Squeal Instability Via Surrogate Modelling[END_REF][START_REF] Renault | Experimental Investigations for Uncertainty Quantification in Brake Squeal Analysis[END_REF][START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF]. These studies remain complex as the different types of uncertainties must be identified and characterized, before using advanced and costly tools to take them into account.

Meta-modeling is now a classical method to propagate uncertainty for an expensive numerical model as it requires a small number of evaluations for its construction. From a few simulations, a mathematical model that surrogates the real model is constructed. Different methods exist to do so, chosen depending on the description of the uncertainty. For random parameters, the polynomial chaos expansion (PCE) [START_REF] Sarrouy | Piecewise Polynomial Chaos Expansion With An Application to Brake Squeal of a Linear Brake System[END_REF][START_REF] Nechak | Robustness of Stochastic Expansions for the Stability of Uncertain Nonlinear Dynamical Systems-Application to Brake Squeal[END_REF], the perturbation method [START_REF] Nobari | Statistics of Complex Eigenvalues in Friction-Induced Vibration[END_REF], or the FORM method [START_REF] Lü | Stability Analysis and Improvement of Uncertain Disk Brake Systems with Random and Interval Parameters for Squeal Reduction[END_REF] can be employed. The fuzzy logic may also be considered [START_REF] Do | A Global Strategy for the Stability Analysis of Friction Induced Vibration Problem With Parameter Variations[END_REF][START_REF] Lü | Uncertainty Quantification of Squeal Instability Under Two Fuzzy-Interval Cases[END_REF]. An other approach consists in the creation of a surface response function with regression [START_REF] Lü | Brake Squeal Reduction of Vehicle Disc Brake System With Interval Parameters by Uncertain Optimization[END_REF][START_REF] Treimer | Uncertainty Quantification Applied to the Mode Coupling Phenomenon[END_REF], the kriging method [START_REF] Nobari | Uncertainty Quantification of Squeal Instability Via Surrogate Modelling[END_REF][START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF][START_REF] Nechak | Sensitivity Analysis and Kriging Based Models for Robust Stability Analysis of Brake Systems[END_REF], support vector machine, neural networks, etc.

However, different types of uncertainties might be present in the same system and the need for hybrid methods that allow different descriptions of the uncertainty in the same surrogate model emerges. Hence, Lu and Yu associate the interval approach with the response surface method (RSM) [START_REF] Lü | Optimization Design of a Disc Brake System With Hybrid Uncertainties[END_REF] for optimizing a brake system regarding squeal. In Ref. [START_REF] Lü | Uncertainty Quantification of Squeal Instability Under Two Fuzzy-Interval Cases[END_REF], they associate the RSM method with the fuzzy logic. In both cases, a RSM model of one eigenvalue of the system is created. This model depends on the interval (resp. fuzzy) variables, and then, the uncertainty is propagated through the RSM. In Ref. [START_REF] Denimal | A Novel Hybrid Surrogate Model and Its Application on a Mechanical System Subjected to Friction-Induced Vibration[END_REF], Denimal et al. presented a new hybrid surrogate model associating kriging with PCE. It enables to deal with both random and parametric variables. The method was tested and validated on a 4-dof phenomenological model to predict its stability. If promising results were obtained, some limitations were raised. First, a high number of learning points were required to build the hybrid meta-model (about 3000 for three parameters). Thus, its application on a full finite element model (FEM) remains an open question as the number of points is limited due to large simulation time. Additionally, when dealing with real FEM of industrial brake systems, one observes a high modal density together with numerous unstable modes and complex coupling phenomena. These make the modal identification complex, and modes can only be tracked when unstable [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF]. This is a major limitation as it means the functions (here the eigenvalues) to be surrogated are only partially defined and sometimes on a small portion of the design space. More concretely, on the total number of available points, only a part of them can be used to create the surrogate model of an eigenvalue, reducing drastically the number of training points. It raises the questions of the ability of the hybrid surrogate model to deal with partially defined problems with a low number of points in the experimental design. So, this work can be seen as an extension of Ref. [START_REF] Denimal | A Novel Hybrid Surrogate Model and Its Application on a Mechanical System Subjected to Friction-Induced Vibration[END_REF] not only to demonstrate the ability and validity of the hybrid surrogate model to predict the stability of industrial FEM of brake systems but also to discuss its limitations.

The automotive brake system is subjected to three uncertain parameters, namely the friction coefficient at the pad/disc interface and two mass additions on the caliper. The first one varies from one braking event to an other and is largely related to the environmental conditions. For this reason, it is commonly described by a probability density function (PDF) [START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part I: Mechanics of Contact and Friction[END_REF][START_REF] Soom | Interactions Between Dynamic Normal and Frictional Forces During Unlubricated Sliding[END_REF][START_REF] Oberst | Statistical Analysis of Brake Squeal Noise[END_REF]. Different PDFs will be considered to demonstrate its high influence and the necessity to identify it correctly. A uniform distribution but also beta distributions will be considered as they are more likely to fit a large variety of experimentally measured distributions.

The two mass additions on the caliper correspond to what one can call design parameters and hence are modeled with an interval. From a practical point of view, in a design process one will look for the best choice of mass addition considering the uncertainty on the friction coefficient. The hybrid surrogate method associating the PCE and the kriging is adopted to propagate both uncertainties simultaneously and to predict each unstable eigenvalue of the system. From the hybrid formulation, a direct access to the stochastic moments of each eigenvalue with regard to the design parameters is assessed without any additional computation. Different probability laws are considered for the friction coefficient, and in each case, the hybrid meta-modeling method is able to predict correctly the evolution of each unstable mode. The results clearly show the strong sensitivity of the brake system stability to the friction coefficient probability law.

Mechanical System Under Study

In this section, the mechanical system as well as the uncertain parameters are presented. Then, the numerical strategy used to predict the squeal propensity is briefly reminded.

Finite Element

Model and Uncertain Parameters. The mechanical system studied is a FEM of an industrial automotive brake system with a floating caliper technology. It is represented in Fig. 1(a). During a breaking action, a hydraulic pressure makes move the piston which pushes the inner pad against the disc. The reaction forces bring back the caliper and the outer pad, which also tackles against the disc. The contact between the pads and the disc creates friction. The contact is modeled with a linear penalty enforcement law and with a Coulomb's law for the friction. For the sake of consistency, the implementation of the FEM and the influence of the internal contacts as well as its validation are not presented here, but the interested reader is invited to refer to Refs. [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF][START_REF] Denimal | Squeal Analysis Based on the Effect and Determination of the Most Influential Contacts Between the Different Components of an Automotive Brake System[END_REF] for more details.

The friction coefficient μ at the pad/disc interface has a high influence on the stability of the brake. Moreover, this parameter is often not well known and varies from one braking action to another one and largely depends on the environment (temperature, humidity, surface, etc.) [START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part I: Mechanics of Contact and Friction[END_REF][START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part II: Dynamics and Modeling[END_REF][START_REF] Butlin | Friction-Induced Vibration: Quantifying Sensitivity and Uncertainty[END_REF][START_REF] Tison | Improvement in the Predictivity of Squeal Simulations: Uncertainty and Robustness[END_REF][START_REF] Kinkaid | Automotive Disc Brake Squeal[END_REF][START_REF] Ouyang | Numerical Analysis of Automotive Disc Brake Squeal: A Review[END_REF]. For this reason, the latter is often considered as variable to indirectly take into consideration these environmental variabilities. Moreover, results shown in Refs. [START_REF] Soom | Interactions Between Dynamic Normal and Frictional Forces During Unlubricated Sliding[END_REF][START_REF] Oberst | Statistical Analysis of Brake Squeal Noise[END_REF] demonstrate that the friction coefficient has a statistical distribution and can be described by a probability density function. When a brake exhibits a too high squeal propensity, structural modifications are investigated to improve the brake with regard to the squeal. From an industrial point of view, a classical structural modification consists in the addition of mass at the extremities of the bracket as it is illustrated Fig. 1(b). These modifications have a high impact on the squeal propensity of the brake system [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF] as they directly modify the modes that are involved in the mode coupling instabilities.

In the present study, these three parameters are considered as uncertain and will be modeled with the hybrid surrogate model. First of all, experimental data on the distribution or evolution of the friction coefficient was not accessible for the example of automotive brake system proposed in our study. This explains why the friction coefficient is taken as random. From the authors' knowledge, the number of available experimental studies in the literature about the distribution law of the friction coefficient in industrial brakes is low. From Ref. [START_REF] Bergman | Influence of Disc Topography on Generation of Brake Squeal[END_REF], the range of variation is taken in [0.3, 0.8].

To illustrate the importance of identifying correctly the friction distribution, three laws are considered here, illustrated in Fig. 2, namely a uniform distribution, which can be seen as "neutral" and two beta laws of shape parameters (2,8) and (8,2) are considered. The choice of using beta distributions is motivated by the idea that they are more likely to fit experimental distribution than a Gaussian one. Indeed, in Refs. [START_REF] Ibrahim | Friction-Induced Vibration, Chatter, Squeal, and Chaos Part I: Mechanics of Contact and Friction[END_REF][START_REF] Soom | Interactions Between Dynamic Normal and Frictional Forces During Unlubricated Sliding[END_REF][START_REF] Oberst | Statistical Analysis of Brake Squeal Noise[END_REF], the distributions observed for the friction coefficient are non-Gaussian. The shape parameters are chosen arbitrarily in this study for illustration. Once the random parameter is generated, a linear transformation on [0.3, 0.8] is applied to get the corresponding physical friction coefficient. On the other hand, two small masses are added on each part of the bracket and can vary between 0 g and 250 g. They correspond to design parameters, and so they are not described by a probability law. They are modeled in the FEM by a modification of the density of a few elements of the caliper, represented in Fig. 1(b).

Complex Eigenvalue Analysis.

The squeal propensity of a brake system is determined from the complex eigenvalue analysis (CEA). It consists in the analysis of the stability of the non-linear equilibrium position.

First of all, the complete nonlinear equation of the automotive brake system can be written as follows:

M Ẍ + C Ẋ + KX + F nl X ( ) = F p (1)
where M, C, and K are the structural mass, damping and stiffness matrices of the FEM, respectively. F p is the external forces applied on the system (i.e., the pressure on the piston) and F nl is The nonlinear sliding equilibrium position U S of the system is the solution of KU S + F nl (U S ) = F p . To study the stability of the solution, a small perturbation of the equilibrium position is considered and the corresponding eigenvalue problem is:

λ 2 M + λC + (K + J nl ) Ψ = 0 (2)
where J nl is the Jacobian matrix of the non-linear forces F nl at the point U S . The eigenvalues λ j = a j + iω j are complex and ω j corresponds to the pulsation of the mode. The mode shape is given by Ψ j . The equilibrium position is considered as stable if all the real parts are negative and is considered as unstable if at least one real part is positive.

Hybrid Surrogate Method

In this section, the hybrid method used is presented. The objective of the surrogate model is to predict the evolution of the eigenvalues with respect to the different uncertain parameters. The influence of the friction coefficient, which is described by a PDF, is modeled with the PCE, whereas the two masses are taken in consideration with a kriging surrogate model.

Mathematical Descriptions

Polynomial Chaos Expansion.

Let's consider that an eigenvalue Λ is a random variable that depends on the random friction coefficient ξ μ . Then, the Generalized Polynomial Chaos (GPC) allows to approximate Λ with a polynomial expansion where the polynomials are orthogonal with respect to the PDF of ξ μ [START_REF] Wiener | The Homogeneous Chaos[END_REF][START_REF] Sudret | Global Sensitivity Analysis Using Polynomial Chaos Expansions[END_REF]. These polynomials are given by the Askey scheme [START_REF] Xiu | Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos[END_REF]. In the case of a uniform law, Legendre polynomials are employed whereas for the case of a Beta distribution, Jacobi polynomials are used. This series is convergent to L 2 sense. For numerical reasons, this series is truncated at the order P. It writes:

Λ(ξ μ ) = P k=0 α k Φ k (ξ μ ) (3)
where (Φ k ) is the polynomial basis and (α k ) is the unknown deterministic coefficient. They can be computed with different nonintrusive techniques such as the regression method [START_REF] Sudret | Global Sensitivity Analysis Using Polynomial Chaos Expansions[END_REF]. Once the coefficients (α k ) are determined, the estimation of the stochastic moments is straightforward. Indeed, the mean of Λ is given by E[Λ] = α 0 and the variance by

σ 2 Λ = P k=0 α 2 k Φ k 2 L 2 .
3.1.2 Kriging. Let's consider a function α that depends on the vector x of variables x M1 and x M2 , the kriging approximates it with the following [START_REF] Forrester | Engineering Design Via Surrogate Modelling: a Practical Guide[END_REF][START_REF] Roustant | Dicekriging, Diceoptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization[END_REF]:

α(x) = K(x) = q i=0 g i (x)β i + Z(x) (4)
where the first part is the regression part and corresponds to the trend of the kriging meta-model. It corresponds to the projection on q functions that are often polynomials of low order. The second part, namely Z(x), is a zero-mean Gaussian process characterized by its covariance Cov

[Z(x), Z(x ′ )] = σ 2 k(x, x ′ )
, where σ 2 is the process variance and k(x, x ′ ) is the correlation function between two input points x and x ′ . k(x, x ′ ) depends on the Euclidean distance h = ‖xx ′ ‖ 2 and on a vector of hyperparameters θ to be optimized.

If the input dimension is superior to 1, the correlation function is constructed with a product of univariate correlation functions.

Hybrid Formulation.

Each eigenvalue of the stability analysis depends on the friction coefficient μ and on the caliper mass modification. The hybrid surrogate model starts from the Eq. ( 3) and considers that the PCE coefficients depend on the mass modification, represented by the variable x. It would be numerically impossible to create a PCE for all the mass modifications; hence, this epistemic uncertainty is modeled with a kriging meta-model. For the m-th eigenvalue, it comes [START_REF] Denimal | A Novel Hybrid Surrogate Model and Its Application on a Mechanical System Subjected to Friction-Induced Vibration[END_REF]:

λ (m) (x, ξ μ ) = P k=0 q i=0 g (m) i (x)β (k),(m) i + Z (k),(m) (x) Φ (m) k (ξ μ ) (5)
To create this meta-model, an experimental design (ED) is generated. An ED corresponds to a set of input and output, it is used to compute the parameters characterizing a meta-model (i.e., the regression coefficient for the PCE, and the regression coefficient as well as the hyper-parameter for the kriging). For the input space, two sets of points are generated. The first one is denoted Ξ = ξ (1) μ , . . . , ξ (N) μ and is composed of N values of the random friction coefficient ξ μ . The second one is denoted A = x (1) , . . . x (M) and is composed of M mass couple values (M 1 , M 2 ). The final ED is obtained by combining these two sets in a full factorial experiment (i.e., N × M configurations). The output set is obtained by evaluating the function to be approximated. This is done by performing a CEA for the N × M system configurations. It is worth emphasizing here that the formulation of the hybrid meta-model makes impossible to generate one ED with the three parameters together as for each x, it is necessary to construct a PCE from N evaluations. The objective is to create two surrogate models for each eigenvalue m: one for the real part, denoted M (m) a and one for the imaginary part, denoted M (m) ω . It is thus necessary to track each eigenvalue over the variations of the system configurations. This is done by using a MAC criterion to track the associated mode shapes [START_REF] Nobari | Uncertainty Quantification of Squeal Instability Via Surrogate Modelling[END_REF][START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF][START_REF] Denimal | A Novel Hybrid Surrogate Model and Its Application on a Mechanical System Subjected to Friction-Induced Vibration[END_REF][START_REF] Denimal | Kriging Surrogate Models for Predicting the Complex Eigenvalues of Mechanical Systems Subjected to Friction-Induced Vibration[END_REF]. However, as demonstrated in Ref. [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF], the MAC criterion is not enough to track the mode at the bifurcation point. Indeed, large variations of the mode shapes are observed at the bifurcation point, making the tracking difficult if the number of points is limited. Moreover, the two complex modes involved in a mode coupling converge to the same real mode at the bifurcation point. It is then impossible, based on a MAC criterion, to differentiate the two modes before the Hopf bifurcation. Hence, the meta-models are created only when the real part of their eigenvalue is non-zero and are valid only on this space. This has a substantial impact on the ED as on the total number of available points, only a few will be available for the surrogate model creation. Moreover, each mode has its own experimental design that corresponds to the part of the total experimental design where the mode is (un)stable (i.e., with a strictly positive or negative real part). Because it does not correspond to a pavement of the input space, a third meta-model is built and acts as a test function to check if the real part is equal to zero or not. The latter is denoted

M (m)
a0 and is built by completing artificially the experimental design of M (m) a with a zero output at the points where the mode was not tracked [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF]. This third meta-model takes zero values when the mode has an eigenvalue with a zero real part, and a non-zero value otherwise (equal to the real part of the eigenvalue), and depending on its value, one can know if the considered mode is stable or not, and so if the surrogate models are in their validity area or not. So, before evaluating a meta-model M (m) a or M (m) ω at a new point (x 0 , ξ 0 ), the test-meta-model M (m) a0 is evaluated to ensure that the two meta-models are in their definition space. If not, one can consider the mode is stable and has a zero-real part. A small threshold equal to 5 is taken in practice to avoid any oscillatory effect. . Once all the PCE are constructed, for each PCE coefficient, a kriging meta-model of the PCE coefficient α (m) j is build based on the available data (i.e., the M values where the mode is unstable). The general workflow is illustrated in Fig. 3.

Application for Squeal Prediction

The section here is devoted to the use of the hybrid surrogate model for the prediction of the squeal propensity of the brake system. In a first time, the construction of the experimental designs is explained, then the meta-models are created and finally the hybrid meta-model is used for the prediction of squeal.

4.1 Construction of the Experimental Designs. As explained previously, three uncertain parameters are considered in the present study: the friction coefficient and two masses added to the bracket. Since the friction coefficient is unknown, it is described by a random law. Three laws are considered here, namely the uniform law and two beta laws of shape parameters (2,8) and (8,2) represented in Fig. 2. The two masses M 1 and M 2 can vary between 0 g and 250 g. A first experimental design composed of 50 couples (M 1 ,M 2 ) is generated with a latin hypercube sampling (LHS) for the kriging, and a second experimental design composed of 15 points is generated by LHS for the friction coefficient for the PCE. Finally, 50 × 15 = 750 input points are used for the hybrid-meta-model creation and the corresponding CEA are Figure 4 shows the stability analysis for the three different laws. Each mode is tracked with an MAC criterion and identified by one color in the complex plans. A large number of unstable modes appear for each configuration with potentially complex crossover phenomena between the different unstable modes. These results show without any ambiguity the complexity of the phenomena involved and consequently the potential difficulty to create an accurate hybrid surrogate model for predicting the stability of a FEM automotive brake system. The influence of the law on the squeal propensity is obvious, and the unstable modes as well as their evolution ranges are also impacted. For example, the mode at 4050 Hz is unstable only when the friction coefficient follows a B(2,8) law. Or, the mode at 4300 Hz has a real part that can reach about 2000 when the friction coefficient follows a uniform law or a B(8,2), but reaches only 1100 for a B(2,8) law. In the following, only the modes with a large enough experimental design are considered for the creation of the hybrid meta-model (i.e., more than 50 points available). This is motivated by two things: first, to have a surrogate model of good quality one needs enough training points; second, if the mode is unstable only on a small part of the experimental design, then one may neglect its unstable behavior in a design process. The threshold of 50 points is arbitrary here, to ensure enough points in the training set. The different modes retained are summarized in Table 1, where for each mode one gives its frequency, its nature, the maximum of the real part (minimum if stable) and the number of points in the experimental design for each law. Again, one can see that the uncertainty related to the friction coefficient has a large impact on the system stability: some modes are unstable only for one PDF, the stability area depends on the PDF (directly related to the number of points in the ED) and variation of the maximum (or minimum) value of the real part of the eigenvalue.

Construction and Validation of the Meta-Models.

Once the experimental designs are generated, the hybrid meta-models are constructed in two steps. First, the PCE is used to estimate the evolution of the frequencies and the real parts of the eigenvalues with respect to the friction coefficient μ. Thus, for each couple (M 1 ,M 2 ), a PCE is constructed for each retained mode. For the uniform law and the beta laws, Legendre and Jacobi polynomials are used, respectively. The degree of each chaos is chosen to be the same for all the modes in one case (same law and same surrogate model) and to obtain the best quality of prediction. To do so, the average relative error between the PCE predictions and the exact values at the points of the ED is minimized. These degrees are summarized in Table 2. For each couple (M 1 ,M 2 ), the number of available points in the experimental design related to the PCE might be small, especially if the mode is mostly stable on the considered interval which explains the low degrees retained for the PCE. As an illustration, the PCE prediction for the real part of the fifth eigenvalue is given in Fig. 5(a). Red points correspond to the training set and the black line to the PCE prediction. One can see the good accuracy of the prediction. Similar results are obtained for the other modes.

Once the PCE are constructed, for each PCE coefficient, a metamodel with the ordinary kriging method is built with a Matérn 5/2 correlation kernel. The kriging predictions of the first two coefficients of the PCE of the M a0 surrogate model of the mode 3 for the three laws are given in Fig. 6. The red points are the learning points, and the surface is the kriging prediction. The influence of the chosen law is quite obvious, the surfaces are complex and vary from one law to another. In the case of the coefficient 0, which is directly related to the mean, the surfaces have similar shapes but the numerical values are different. Indeed, the maximum is equal to 90, 110, and 82 for the uniform law, the B(8,2) law and the B(2,8) law, respectively. For the coefficient 1, both the shapes and the values are strongly influenced by the law. From this, its comes that the mean of the real part of the eigenvalue of the mode at 2850 Hz has a similar evolution (but more or less higher values) for different PDF of the friction coefficient, but that the variance will be strongly impacted. The different metamodels have been validated by comparison of the hybrid metamodel prediction to reference values. To do so, a set of 45 reference points has been generated based on additional simulations. For each mode, the PCE-kriging predictions are compared to these reference values. The number of validation points depends not only on the considered mode but also on the probability law as the surrogate models are not valid on the full design space. As an illustration, the comparison between the predictions for the real part of the eigenvalue of the fifth mode and the reference values are given in Fig. 5(b) where red circles are the reference values and black crosses the predictions. The accuracy of the surrogate model is obvious. This comparison is summarized for each real and imaginary parts of each eigenvalue and for each law in Table 3, where the average relative error e for each mode is given by

e = 1 n n k=1 y (pred) k -y (ref ) k y (ref ) k ( 6
)
where y (pred) k is the prediction of the considered meta-model at the kth validation point and y (ref ) k the reference value. To be noted that the quantities, y (pred) k and y (ref ) k are associated with the real or imaginary parts of the eigenvalue for one specific mode. n corresponds to the number of points in the experimental design for which the real part of the associated mode is superior to zero. Considering more specifically the average relative error e for imaginary parts, all the results indicate unambiguously the relevance and validity of the PCE-kriging predictions. Considering results for real parts, the level of error e is always low giving confidence in the different hybrid surrogate models except for some limited specific cases. These specific cases for which an important error e is detected correspond to modes for which we have both a very small number of points n to calculate e as well as positive real parts very close to zero which has the consequence of increasing the relative error e.

Squeal Prediction.

Once the meta-models are constructed, it is possible to predict the evolution of the different eigenvalues.

The main interest of the present methodology is the estimation of the PCE coefficients with the kriging. Since the PCE coefficients are directly related to the stochastic moments, it is possible from the kriging to predict straightforwardly the evolution of the mean and the variance of the eigenvalues without any MCS as it would have been the case if a unique kriging meta-model was built for the three parameters [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF]. Thus, means and variances with respect to the masses M 1 and M 2 can be directly evaluated. For example, the evolution of the mean and the variance of the real part of the modes 3 and 4 are given in Figs. 7 and8 for the three different laws.

Considering the stability behavior of the mode 3, looking at the average of the real part of the eigenvalue (see in Fig. 7 upper line), the behavior and the evolution is similar over the three different PDF of the friction coefficient. Indeed, the mode is mostly stable, and when M 1 and M 2 increase simultaneously, the mode becomes unstable. This unstable area is slightly smaller for a B(2,8) law. If the evolution is similar, the maximum mean value of the real part is impacted by the friction PDF and larger mean values are expected with a B(8,2) law for example. On the opposite, the behavior of the variance is strongly impacted by the PDF (see in Fig. 8 upper line). For a uniform law a large variance is observed on the whole unstable part, for a B(8,2) law the variance is low and the maximum is observed at the bifurcation area around M 1 = M 2 = 50 g and for a B(2,8) law a large variance is observed only at the most unstable area. The extremely large variance for the uniform law is of importance, as it might imply a change in the stability behavior of the mode depending on the friction coefficient value, and so depending on the environmental conditions a squeal instability at this frequency can be observed or not. But if the friction coefficient is described by a B(2,8) or a B(8,2) law, the variance on the unstable part is low and the mode will always be unstable whatever the friction coefficient value and the environmental conditions. Looking at the stability behavior of the mode 4, the evolution of the mean of the real part is similar but it impacts the stability behavior of the mode (see in Fig. 7 bottom line). First, the average real part is always much higher when M 1 = 50 g. Then, a more or less large stability area is observed for low values of M 2 and large values of M 1 . This stability area is almost non-existent for a uniform friction law, and much larger when the friction coefficient follows a B(2, 8) law. When the two masses have higher values, the mode is also unstable but with lower values of the real part in average for all cases. In this case, the PDF that describes the friction coefficient has an impact on the average stability behavior of the mode. As for mode 3, the variance is strongly impacted by the PDF of the friction law (see in Fig. 8 bottom line). For the uniform law, a large variance is observed in the bottom left part (area where M 1 + M 2 < 250 g). For low values of M 2 and large values of M 1 , it means that the mode is sometimes stable and so that squeal event might take place, or not, depending on the environment. For the B(8,2) case, the variance is always very low and the stability behavior of the mode is not impacted. For the B(2,8) case, the variance is always low except at M 1 = 50 g and M 2 > 50 g. So, if in average the mode has a higher real part at M 1 = M 2 = 50 g than at M 1 = M 2 = 200 g (and could be considered as "more unstable"), it is also more fleeting considering the variance as the variance is extremely low at M 1 = M 2 = 200 g.

Finally, it is possible to access directly the mean number of instability and the sum of the mean of the positive real parts with regards to the masses M 1 and M 2 , they are given in Fig. 9 for each law. The influence of the law for the friction coefficient is obvious. Indeed, when a uniform law is considered, from 2.5 to 5 instabilities are observed in mean, whereas for a B(8,2) law, up to six instabilities can be observed in average. Moreover, the location of the maximum and minimum number of instabilities is different depending on the law. Indeed, for a uniform law, a mass M 1 null and M 2 equal to 250 g gives the lowest number of instabilities, whereas for a B(2,8) law, the minimum is reached for both M 1 null and M 2 inferior to 80 g. Also, the shape of the limits between the average number of instabilities (see black lines in Fig. 9) is completely different from one case to another.

A similar analysis can be done for the sum of the mean of the positive real parts, given in Figs. 9(d), 9(e), and 9( f ). Indeed, the mean values are highly impacted by the distribution: see the different ranges of variations in Figs. 9(b), 9(d), and 9( f ). For a B(2,8) law, the highest value is equal to 700, whereas for a B(8,2), it is equal to 2600. The locations of the maximum and minimum are also impacted. Indeed, for a B(2,8) law, the minimum is reached for M 1 = 50 g and M 2 = 40 g, whereas for a uniform law, it is reached for M 1 = 75 g and M 2 = 240 g.

Results obtained here demonstrate the complex influence of the design parameters on the system stability which illustrates the necessity to conduct deep studies to end in a correct understanding of squeal phenomenon. The evolution of the different modes is not trivial which complicates the identification of the best set of design parameters for an engineer. Moreover, the identification of the probability laws of the random parameters must be done beforehand and precisely as it has a high impact on the system stability. 4.4 Discussion. The results presented give a concrete use of a hybrid surrogate model for the prediction of stability of an industrial brake system, but also raises some limitations:

size of the ED (here 750): it is quite high and is definitely a limitation related to the use of this method, but this high number is related to different factors. The first reason is the full factorial construction of the ED imposed by the formulation of the hybrid surrogate model as the uncertainties are propagated one after another, and work should be conducted to lift this limitation. However, this high number of points is mostly due to the physical problem under consideration. Indeed, due to the high modal density and the important variation of the mode shapes over the design space, the ED must be dense enough to ensure a robust mode tracking, which implies a minimal number of points. Moreover, as modes can be tracked only when the eigenvalues real parts are non-zero, the ED must be dense enough to ensure that each instability is indeed detected, but it also implies that on the total number of points, only a limited number is finally available to create a surrogate model for a mode. So, when 750 simulations are run, in some cases only 80 points are really available comparison with a kriging with all the parameters: one could naturally consider using one kriging meta-model for the three parameters [START_REF] Denimal | Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method[END_REF] and then perform MCS for the propagation of the random law related to the friction coefficient. Using this strategy represents a reduction of the initial number of CEA to perform. However, the MCS to be performed at the end takes time, several minutes as many modes are present and as the validity of the surrogate models must be tested each time. By using the hybrid surrogate model, this simulation time is removed as the stochastic moments are directly obtained from the PCE coefficients. In some cases, it might be more convenient to have a more important simulation time initially and get directly statistical moments at the end, than having to run again MCS after the surrogate model creation. One last advantage, not illustrated here as only one parameter is random, is that the Sobol indices can also be directly accessed from the PCE coefficients. computational time: except the simulations related to the CEA, the most computationally expensive part of the method is the mode tracking part, with the computation of many MAC matrices over the ED and represents a few minutes in the present case. The construction of the hybrid meta-models represents only 1.5 s on average, showing the efficiency of the method. Once the surrogate models are constructed, obtaining the evolution of the mean and variance on a 50 × 50 grid (Fig. 7) takes only 0.13 s for a mode, when a few minutes would be required with a MCS on a single kriging surrogate model.

Conclusion

The present study proposes the application of a hybrid surrogate method for the prediction of the stability behavior of an industrial brake system subjected to many instabilities. The main objective is to demonstrate the validity and the interest of hybrid methods when a system is put through different types of uncertainties. Indeed, for many application epistemic and aleatory uncertainties are present, and modeling them at the same time in the same surrogate model is not necessarily the best representation.

The CEA methodology is employed for the prediction of the stability behavior of the finite element model of the brake system and three uncertain parameters are considered. First, the friction coefficient is described by three different probability laws, and two design parameters corresponding to two small masses that might be added to the caliper are described by an interval of possibility. The friction coefficient is modeled via a PCE whereas the two design parameters are modeled through the kriging method. Once the meta-models are constructed and validated, the analysis of the stability of the brake system can be performed. The work demonstrates the applicability of such an approach for large models where computational times are important and where the construction of the learning sets is limited due to the mode tracking, which limits drastically the number of available points to construct the meta-models.

From the PCE coefficients, it is possible to access without any additional computations (i.e., no MCS) to the evolution of the mean and the variance of each eigenvalue regarding the design parameters allowing for a stochastic description of the system. This would not be possible if a unique surrogate model for all the uncertain parameters was created, since MCS simulations would be performed to propagate the uncertainty related to the random parameters. The approach used here may require more simulations to create the surrogate models, but no cost is added then for the UQ part. Considering the global process, this strategy might be of interest in some cases when stochastic properties are required.

Moreover, the strong impact of the probability law chosen for the friction coefficient on the system stability is also observed and demonstrates the necessity to identify it in advance.

Fig. 1 (

 1 Fig. 1 (a) Exploded view of the brake and (b) mass modification on the caliper

Fig. 2

 2 Fig. 2 Considered probability laws for the random parameter associated with the friction coefficient μ: (a) uniform, (b) B(8,2), and (c) B(2,8)

Finally, for

  each mode m and for each mass couple (i.e., for each x (k) ∈ A), a PCE is constructed based on the N model evaluations for which the mode is unstable. The PCE coefficients are the α (k),(m) j j∈ 0,P { }

Fig. 3

 3 Fig. 3 General workflow of the PCE-Kriging construction

Fig. 4

 4 Fig. 4 Experimental design from the mode tracking for the different laws: (a) uniform law, (b) B(8,2) law, and (c) B(2,8) law

Fig. 5

 5 Fig. 5 Real part of the eigenvalue of the mode 5: (a) comparison of the PCE prediction (black line) to reference points (red dot) for (M 1 , M 2 ) = (212.1, 137.3) g and (b) comparison of the PCE-Kriging prediction (black cross) to different validation points (red circle) for different values of μ, M 1 , and M 2 (Color version online.) 021006-6 / Vol. 144, APRIL 2022 Transactions of the ASME

Fig. 6

 6 Fig. 6 Kriging prediction of the first two PCE coefficients for the real part of the mode 3 for the three laws with respect to the masses M 1 and M 2 -ED (•): (a,d) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law (Color version online.)

Fig. 7

 7 Fig. 7 Evolution of the mean of the real part with regard to the masses M 1 and M 2 for the modes 3 (top) and 4 (bottom) for the ((a,d)) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law

Table 1

 1 Modes involved in mode coupling for the different experimental designs with their frequency, real part, and the number of points in each experimental design for each law

				Uniform law	B(8,2) law	B(2,8) law
	Mode no.	Frequency (Hz)	Nature	Real part	No. pts ED	Real part	No. pts ED	Real part	No. pts ED
	1	879.7	Unstable	-	-	27.5	112	22.9	100
	2	1855.1	Unstable	248.3	403	248.3	475	217.0	324
	3	2853.2	Unstable	131.1	168	131.1	220	89.5	80
	4	3460.8	Unstable	124.1	601	387.8	696	96.5	459
	5	3863.5	Unstable	670.2	741	774.2	750	406.4	678
	6	4395.7	Unstable	1948.9	739	1948.9	750	1194.9	674
	7	4879.5	Unstable	249.9	409	-	-	82.6	275
	8	5541.1	Unstable	-	-	-	-	54.9	144
	9	879.7	Stable	-	-	-27.5	112	-22.9	100
	10	1855.1	Stable	-248.3	403	-248.3	475	-217.0	324
	11	2853.2	Stable	-131.1	168	-131.1	220	-89.5	80
	12	3460.8	Stable	-124.1	601	-387.8	696	-96.5	459
	13	3863.5	Stable	-670.2	741	-774.2	750	-406.4	678
	14	4395.7	Stable	-1948.9	739	-1948.9	750	-1194.9	674
	15	4879.5	Stable	-249.9	409	-	-	-82.6	275
	16	5541.1	Stable	-	-	-	-	-54.9	144

Table 2

 2 Polynomial chaos order for the different laws and the different surrogate models

	Law	M ω	M a	M a0
	Uniform	4	5	7
	B(8,2)	5	5	3
	B(2,8)	3	3	3

Table 3

 3 Average relative error e of the PCE-kriging predictions at validation points on the eigenvalues estimation for the different laws

		Error real part		Error imaginary part
		Uniform B(8,2)	B(2,8) Uniform B(8,2)	B(2,8)
	Mode 1	-	8.51E-2 1.04	-	2.57E-4 5.07E-4
	Mode 2	4.45E-2 1.95E-1 2.96E-1 1.01E-3 9.53E-3 2.34E-3
	Mode 3	7.35E-2 0.42E-2 2.97E-1 6.38E-4 2.52E-4 2.1E-3
	Mode 4	2.63E-1 1.31E-1 4.82E-1 1.12E-3 8.12E-4 7.00E-4
	Mode 5	1.42E-2 2.59E-2 3.75E-2 1.74E-4 7.82E-4 4.17E-4
	Mode 6	7.65E-2 6.28E-3 1.11E-1 1.32E-3 1.14E-3 1.33E-3
	Mode 7	4.88E-1	-	7.65E-2 4.10E-4	-	8.28E-5
	Mode 8	-	-	4.48E-1	-	-	2.87E-4
	Mode 9	-	8.51E-2 1.04	-	2.57E-4 5.07E-4
	Mode 10 4.51E-2 2.55E-1 2.96E-1 1.04E-3 9.53E-3 2.39E-3
	Mode 11 7.35E-2	-	2.97E-1 6.38E-4	-	2.14E-3
	Mode 12 3.61E-1 2.73E-1 4.69E-1 5.36E-2 1.66E-3 7.26E-4
	Mode 13 1.42E-3 2.59E-2 3.44E-2 1.74E-4 7.82E-4 4.34E-4
	Mode 14 7.65E-2 6.28E-3 1.1E-1	1.32E-3 1.15E-3 1.33E-3
	Mode 15 4.88E-1	-	7.65E-2 4.10E-4	-	8.28E-5
	Mode 16	-	-	4.47E-1	-	-	2.87E-4

Journal of Vibration and Acoustics

APRIL 2022, Vol. 144 / 021006-7

Acknowledgment

This work was achieved within Stellantis Stellab program-OpenLab Vibro-Acoustic-Tribology@Lyon (VAT@Lyon). J.-J. Sinou acknowledges the support of the Institut Universitaire de France. Fig. 9 Mean number of unstable modes (top) and sum of the mean of the positive real parts (bottom) with regard to the masses M 1 and M 2 for (a,d) the uniform law, (b,e) the law B(8,2), and (c,f) the law B(2,8). Black lines in (a,b,c) denote the contours of integer numbers of instabilities. Journal of Vibration and Acoustics APRIL 2022, Vol. 144 / 021006-9

Conflict of Interest

There are no conflicts of interest.


Ma -Olivier Titaud -Arthur

Vidard -François-Xavier Le Dimet

Jianwei Ma

Olivier Titaud
email: olivier.titaud@imag.fr

Arthur Vidard
email: arthur.vidard@imag.fr

François-Xavier Le Dimet
email: francois-xavier.le-dimet@imag.fr

Spatio-temporal structure extraction and denoising of geophysical fluid image sequences using 3D curvelet transforms

Keywords: 3D Curvelets, feature extraction, 3D total variation, data assimilation, video/image sequences, remote sensing

Since several decades many satellites have been launched for the observation of the Earth for a better knowledge of the atmosphere and of the ocean. The sequences of images that such satellites provide show the evolution of some large scale structures such as vortices and fronts. It is obvious that the dynamic of these structures may have a strong predictive potential. Extracting these structures and tracking their evolution automatically is then essential for future forecast systems. In this paper we consider extraction of spatio-temporal geometric structures from image sequences of geophysical fluid flow using threedimensional (3D) curvelet transform and total variation minimization. Numerical experiments on simulated geophysical fluids and real video data by remote sensing show good performances of the proposed method in terms of denoising and edge structural extraction. This work is partially motivated by a sequent application to image sequence assimilation of geophysical fluids.

Extraction de structures spatio-temporelle et débruitage de séquences d'image de fluides géophysiques à l'aide de transformées de courbelettes 3D

Résumé : Depuis quelques décénies de nombreux satellites d'observation de la terre ont été lancésafin d'améliorer nos connaissances de l'atmosphère et de l'océan. Les séquences d'image fournies par de tels satellites montrent l'évolution de structures grandes échelles telles que les syclones et les fronts. Il est évident que la dynamique contenue dans ces structures peuvent avoir un fort potentiel prédictif. Extraire ces structures et suivre leur évolution de façon automatique est donc essentiel pour les futurs systèmes de prévision. Dans cet article on considère l'extraction de structures géometriques spatio-temporelles dans des séquences d'images de fluides géophysiques en utilisant des transformées en coubelettes 3D une minimisation de la variation totale. Des expériences numériques sur des images de fluides géophysiques simulées et des données de vidéo réelles montrent la bonne performance de la méthode proposée en terme de débruitage et extraction de structure.

Mots-clés : Courbelettes 3D , extraction de caractéristiques, Total Variation 3D, assimilation de données, video/séquences d'image, télédetection 1 Introduction

Background

Current weather and ocean forecast systems assimilate [START_REF] Dimet | Variational algorithms for analysis and assimilatoin of meteorological observations: theoretical aspects[END_REF][START_REF] Blum | Data Assimilation for Geophysical Fuilds, Handbook of Numerical Analysis[END_REF] in-situ or remote measurements of the state variables of the modeled flow such as winds, pressure, temperature, salinity,. . . in order to estimate initial conditions. Since several decades many satellites have been launched for the observation of the Earth for a better knowledge of the atmosphere and of the ocean. Geostationary or polar orbiting Earth observation satellites provide image sequences which clearly show the evolution of large scale features such as fronts, edges shapes of vortices and their trajectory (see Fig. 1 and2).

Unfortunately, this kind of information are not yet taken into account in a quantitative way by the operational forecast centers. However the evolution of these structures have a strong predictive potential because they may give some information about the dynamic of the underlying physical system.

Extracting and tracking fluid structures from geophysical image sequences could prove very valuable in the framework of Image Sequences Assimilation (ISA). Indeed, the later aim at improving forecast systems by taking into account dynamical informations contained in the structure evolution that satellite images provide [START_REF] Dimet | Assimilation of Images In Geophysical Models[END_REF][START_REF] Herlin | Coupling models and data: which possibilities for remotely-sensed images?[END_REF][START_REF] Dimet | Assimilation of Images in Numerical Models in Geophysics[END_REF].

In the prior work [START_REF] Ma | Curvelet-based snake for multiscale detection and tracking of geophysical fluids[END_REF], we applied curvelet transform [START_REF] Candès | Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF][START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise singularities[END_REF] to the detection and tracking of characteristic deformable structures in some geophysical fluid

Motivation

Although applications of wavelets to extraction and analysis of geophysical fluid images have become increasingly popular (see e.g. [START_REF] Mallat | Singularity detection and processing with wavelets[END_REF][START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF][START_REF] Farge | Coherent vortex extraction in 3D turbulent flows using orthogonal wavelets[END_REF][START_REF] Turiel | Wavelet filtering to extract coherent vortices from altimetric data[END_REF][START_REF] Yano | Wavelet analysis of simulated tropical convective cloud systems. Part I: Basic analysis[END_REF][START_REF] Chapa | Application of wavelet transform to Meteosatderived cold cloud index data over South America[END_REF] among lots of literature), traditional wavelets ignore the geometric properties of structures and do not exploit the regularity of edges : they are efficient to represent point singularities only. Therefore, wavelet-based structure extraction and denoising becomes computationally inefficient for geometric features with 2D line and 3D surface singularities. Moreover discrete wavelet thresholding could lead to oscillations along edges.

The curvelet transform was proposed by Candès et al. [START_REF] Candès | Ridgelets: A key to higher-dimensional intermittency ?[END_REF][START_REF] Candès | Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF][START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise singularities[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF]. It is a new geometric wavelet transform. The 2D curvelet transform allows an optimal sparse representation of objects with C 2 -singularities. In comparison with wavelets, curvelets show better performances for the analysis of linesingularities. Thus curvelet seems to be more adapted to extract edges and then structures. For a smooth object f with discontinuities along C 2 -continuous curves, the best m-term approximation fm by curvelet thresholding obeys the inequality f -fm 2 2 ≤ Cm -2 (log m) 3 , while for wavelets the decay rate is only m -1 . Surprising performance has been shown in fields of image processing, see e.g. [START_REF] Starck | The curvelet transform for image denoising[END_REF][START_REF] Starck | Astronomical image representation by the curvelet transform[END_REF][START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF][START_REF] Ma | Curvelets for surface characterization[END_REF][START_REF] Ma | Combined complex ridgelet shrinkage and total variation minimization[END_REF]. Recently, the 3D curvelet transform was presented by Ying et al. [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF]. Unlike the isotropic elements of wavelets, the needle-shaped elements of this transform possess very high directional sensitivity and anisotropy (see Fig. 1 for 2D case). Such an element is very efficient for representing vortex edges of geophysical fluids images.

In a previous work [START_REF] Ma | Curvelet-based snake for multiscale detection and tracking of geophysical fluids[END_REF], 2D curvelet transforms have been applied to the detection and tracking of some characteristic deformable structures of geophysical fluid in sea surface temperature satellite image sequences (from NOAA Advanced Very High Resolution Radiometerimage). However such detection do not take into account temporal consistency between successive frames. In this paper, we apply a 3D curvelet transform to fluid flow image sequences in or-INRIA der to extract the evolution of spatio-temporal structures. Then a new curvelet reconstruction by combining the 3D curvelet transform with total variation minimization [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF][START_REF] Ma | Curvelets for surface characterization[END_REF][START_REF] Ma | Combined complex ridgelet shrinkage and total variation minimization[END_REF] is presented to further suppress the artifacts resulting from classical thresholding.

The paper is organized as follow : part II is devoted to the description of 3D curvelet transform and some classical threshold operations; part III describes the TV minimization technique and its combination with 3D curvelet transform: this forms the original part of this paper. Some numerical experiments showing applications of the later method to geophysical fluid flow image sequences are presented in part IV.

3D curvelet transform and thresholding 2.1 2D curvelet transform

In 1999, an anisotropic geometric wavelet called ridgelet was proposed by E. Candès and D. Donoho [START_REF] Candès | Ridgelets: A key to higher-dimensional intermittency ?[END_REF]. It is optimal at representing straight-line singularities. However its application is limited to objects with global straight-line singularities. In order to take into account images with local line or curve singularities, the so-called first generation of Curvelet transform [START_REF] Candès | Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF] was developed using ridgelet transforms applied on sub-partitions of the image. A second generation of curvelet transform based on a frequency partition technique was introduced by the same authors in 2004 [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise singularities[END_REF]. A 2D discrete curvelet transform was implemented in 2004 [START_REF] Candès | Fast discrete curvelet transforms[END_REF]: it is freely available at http://www.curvelet.org.

The 2D continuous curvelet transform can be defined as follow. Let W : [1/2, 2] → R and V : [-1, 1] → R be respectively a radial and an angular windows -in the spatial domain -, supported on [1/2, 2] and [-1, 1] respectively and obeying the admissibility condition

-∞ l=-∞ V 2 (t -l) = 1, t ∈ R, (1)
-∞ l=-∞ W 2 (2 -j r) = 1, r > 0. (2)
For instance, the following Meyer windows [START_REF] Ma | Curvelets for surface characterization[END_REF] satisfy the above conditions:

V (t) =    1 |t| ≤ 1/3, cos[ π 2 ν(3|t| -1)] 1/3 ≤ |t| ≤ 2/3, 0 else. (3)
W (r) = min cos( π 2 ν(5 -6r)), cos( π 2 ν(3r -4)) ( 4
)
where ν is a smooth function satisfying for all x ∈ R :

ν(x) = 0 x ≤ 0, 1 x ≥ 1, and ν(x) + ν(1 -x) = 1. (5) ν(x) =    0 x ≤ 0, 1 2 (sin[ π 2 (2x -1)] + 1) x ≤ 1, 1 else. ( 6
)
RR n°6683 As in wavelet theory, we also introduce a low-pass window W j0 , j 0 ≥ 0 which satisfies : for all r ∈ R

W j0 (r) 2 + j>j0 W 2 (2 -j r) = 1. (7)
Finally, let us denote by f the Fourier transform of an element f of L 2 (R) : it is a function of the frequency variables ξ = (ξ 1 , ξ 2 ) with polar coordinates (r, θ). Curvelets are indexed by three integers j, k and l which respectively stand for scale, location and orientation. Curvelets at the coarsest scale j = j 0 are defined through their Fourier transform by

ϕ j0,0,k (x) = ϕ j0 (x -2 -j0 k), k = Z 2 , (8)
where

φj0 (ξ) = 2 -j0 W j0 (2 -j0 |ξ|). (9)
At this scale, curvelets are then isotropic.

Curvelets at the finest scales. For each scale j > j 0 , we define a window U j in the frequency domain for all ξ = (ξ 1 , ξ 2 ) with polar coordinates (r, θ) by

U j (ξ) = 2 -3j/4 W (2 -j r)V 2 ⌊j/2 ⌋ θ 2π , (10)
where ⌊•⌋ denotes the integer part operator. The support of U j is a polar 'wedge' determined by supp(r 4 shows the first window U 0 in the frequency domain when the radial and angular windows are defined by ( 3) and ( 4) respectively and when ν given by ( 6) [START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF]. We consider a uniform discretisation of the angle coordinate θ ∈ [0, 2π[

→ W (2 -j r)) = [2 j-1 , 2 j+1 ] and supp(θ → V (2 ⌊j/2⌋ θ)) = [-2 -⌊j/2⌋ , 2 -⌊j/2⌋ ]. Figure
θ l = 2π × l × 2 -⌊j/2⌋ , 0 ≤ l < 2 ⌊j/2⌋ , (11)
and the rotation by θ radians matrix

R θ = cos θ sin θ -sin θ cos θ . ( 12
) INRIA ∼ 2 j/2 ∼ 2 j
Figure 5: Frequency plane partition. The shaded area shows an example of a "wedge"on which a curvelet is supported in the frequency plane. At scale j, this wedge has a width of order 2 j and a length of order 2 j/2 .

Note that the number of discrete orientations θ l depends on the scale index j : this number is doubled when scale change from an odd index 2j + 1 to the larger next even one 2j + 2. It is unchanged when scale index is changing from 2j to 2j + 1 (see Fig. 5). At the finest scales j > j 0 , curvelets are defined for

all x ∈ R 2 by ϕ j,k,l (x) = ϕ j R θ (x -x (j,l) k ) ( 13
)
where the waveform ϕ j is the inverse Fourier transform of the frequency window U j and the translation vector x

(j,l) k , k = (k 1 , k 2 ) ∈ Z 2 is defined by x (j,l) k = R -1 θ l (2 -j k 1 , 2 -j/2 k 2 ). ( 14
)
It can be proved [START_REF] Candès | New tight frames of curvelets and optimal representations of objects with piecewise singularities[END_REF] that the above family of curvelets defines a tight frame (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] for a detailed definition of this term): for all f ∈ L 2 (R 2 ) we have (Parseval equality)

j,k,l | f, ϕ j,k,l | 2 = f 2 L 2 (R 2 ) , (15)
and the reconstruction formula

f = j,k,l ϕ j,k,l , f ϕ j,k,l , (16)
where

•, • is the usual inner product in L 2 (R 2 ): f, g = R 2 f g. Remarks:
• Curvelet coefficients in ( 16) are determined by

c j,k,l := ϕ j,k,l , f = 1 (2π) 2 f (ξ)U j (R θ l ξ)e i x (j,l) k ,ξ dξ (17)
• In practice we use Discrete Curvelet Transform (DCT) which is designed for f defined on a Cartesian grid

{(n 1 , n 2 ) ∈ N, 0 ≤ n 1 , n 2 < n ∈ N}. Digital coefficients are defined c D j,k,l = n1,n2 f (n 1 , n 2 )ϕ D j,k,l (n 1 , n 2 ) ( 18
)
where the construction of the digital curvelet waveforms ϕ D is an adaptation of the above continuous transform for Cartesian grids. This lead to a non uniform discretisation of the orientations. For more details and C++ implementation see below and [START_REF] Candès | Fast discrete curvelet transforms[END_REF].

3D curvelet transform

The 3D curvelet transform was proposed by Ying et al. [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF]. A fast application of the 3D curvelet transform to 3D turbulent flows has been presented by Ma et al. [START_REF] Ma | Three-dimensional curvelets for coherent vortex analysis of turbulence[END_REF]. The so-called second-generation curvelet discrete transform used in this paper is based on a frequency partition technique. Similar to 2D problems, we define V (t) and W (r) to be a pair of smooth, non-negative real-valued window functions, which are called angular window and radial window, respectively. V is supported on [-1, 1] and W on [ 1 2 , 2]. The windows should satisfy admissibility conditions described in [START_REF] Candès | Fast discrete curvelet transforms[END_REF]. Without loss of generality, we introduce the lowpass window W 0 for the coarsest scale, which satisfies the condition

W 0 (r) 2 + j>0 W (2 -j r) 2 = 1. ( 19
)
For example, these conditions are satisfied by taking the scaled Meyer windows ( 3) and ( 4) [START_REF] Ma | Curvelets for surface characterization[END_REF]. For each j > 0, the radial window W (2 -j r) smoothly extracts the frequency content inside the dyadic corona 2 j-1 ≤ r ≤ 2 j+1 . The angular windows partition of R 3 into trapezoidal regions obeying frequency parabolic scaling (thickness ≈ lengh 2 ) [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF]. For j ≥ 0, define the window U j (ξ), ξ = (ξ 1 , ξ 2 , ξ 3 ) ∈ R 3 in frequency domain as

U j (ξ) = 2 -3j/4 W (2 -j |ξ|) V (2 ⌊j/2⌋ θ), ξ ∈ R 3 , (20)
where (|ξ|, θ) denotes the polar coordinates corresponding to ξ. The support of U j is a 3D polar wedge shape, i.e., a half circular cone (see Fig. 3 in [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF]).

The system of curvelets is now indexed by three parameters (j, l, k) where j denotes scale, l denotes orientation, and k = (k 1 , k 2 , k 3 ) denotes spatial location. Define the curvelets as

ϕ j,l,k (x) := ϕ j (R θ j,l (x -x (j,l) k )), x = (x 1 , x 2 , x 3 ) ∈ R 3 (21)
where ϕ j (ξ) := U j (ξ), i.e., U j is the Fourier transform of ϕ j , R θ j,l denotes the rotation matrix with angle θ j,l . Figure 6 shows an element of 3D curvelets. Observe that in the spatial domain, ϕ j,l,k is of plate-like shape, which rapidly decays away from a 2 -j by 2 -j/2 cross-section rectangle with center x (j,l) k and orientation θ j,l with respect to the vertical axis in x. The element is smooth within the plate but exhibits oscillating decay in the normal direction of the plate. It obeys a parabolic scaling law between the thickness and length (thickness ≈ lengh 2 ) and directional sensitivity (Orientations = 1/ √ scale). Let µ = (j, l, k) be the collection of the triple index again. The curvelet coefficients are given by In order to have Cartesian coronae, which is based on concentric cubes instead of sphere, the authors of [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF] applied a pseudo-polar grid by a modified window of the form

c µ (f ) : = f, ϕ µ = R 3 f (ξ) φµ (ξ)dξ = R 3 f (ξ) U j (R θ j,l ξ) e i x (j,l) k ,ξ dξ. ( 22
) INRIA (a) (b) (c) (d) (e) (f)
W 0 (ξ) = Φ 0 (ξ), W j (ξ) = Φ 2 j+1 (ξ) -Φ 2 j (ξ), j > 0, (23)
where

Φ j (ξ 1 , ξ 2 , ξ 3 ) = φ(2 -j ξ 1 )•φ(2 -j ξ 2 )•φ(2 -j ξ 3 ), the one-dimensional window φ satisfying 0 ≤ φ ≤ 1, supp φ ⊂ [-2, 2] and φ(r) = 1 for r ∈ [-1/2, 1/2].
As before, φ can be taken to be a scaled Meyer window. The angular window for the l-th wedge can be redefined as

V j,l (ξ) = V (2 j/2 • ξ 2 -α l • ξ 1 ξ 1 ) • V (2 j/2 • ξ 3 -β l • ξ 1 ξ 1 ). (24)
Here (1, α l , β l ) is the direction of the centre line of the wedge. Every Cartesian corona has six components. The windows in the other five components have similar definitions.

Ma et al

Now we define the modified frequency window U as

U 0,0 (ξ) = W 0 (ξ), j = 0, (25)
U j,l (ξ) = W j (ξ) • V j,l (ξ), 0 < j < j f , (26)
and at the finest scale j f , the wavelet-like isotropic element is defined by the frequency window U j f ,0 (ξ) = W j f (ξ). It is clear that U j,l (0 < j < j f ) isolates frequencies near the wedge

(ξ 1 , ξ 2 , ξ 3 ) : 2 j-1 ≤ ξ 1 ≤ 2 j+1 , -2 -j/2 ≤ ξ 2 /ξ 1 -α l ≤ 2 -j/2 , -2 -j/2 ≤ ξ 3 /ξ 1 -β l ≤ 2 -j/2 . ( 27
)
Let L p,j,l (p = 1, 2, 3) be three positive integers satisfying (i) One can not find ξ and ξ ′ such that ξ p -ξ ′ p are multiples of L p,j,l ;

(ii) The volume Λ j,l = L 1,j,l • L 2,j,l • L 3,j,l is minimal.

The two previous conditions guarantee that the data does not overlap with itself during the wrapping process below. Obviously, the U j,l is supported now in a 3D rectangular box of integer size L 1,j,l × L 2,j,l × L 3,j,l .

The discrete curvelets are given by their Fourier formation

φD µ (ξ) = U j,l (ξ) • exp -2πi p=1,2,3 k p ξ p L p,j,l / Λ j,l (28)
for 0 < k p < L p,j,l , p = 1, 2, 3. Analogously, the transform at the coarsest level is defined as

φD 0,0,k (ξ) = U 0,0 (ξ) • exp -2πi p=1,2,3 k p ξ p L p,0 / Λ 0 , (29)
and a similar formula can be obtained at the finest scale by replacing the scale 0 with j f and setting L p,j f = n. Now we can find the Cartesian counterpart of the coefficients in [START_REF] Korotaev | Retrieving ocean surface current by 4-D variational assimilation of sea surface temperature images[END_REF] by

cD µ (f ) = f, ϕ D µ = R 3 W( U j,l (ξ) f (ξ)) e i kj ,ξ dξ. ( 30
)
The algorithm of three-dimensional discrete curvelet transform is summarized as follow:

1. Apply the 3D FFT and obtain Fourier samples f (ξ), -n/2 ≤ ξ < n/2, ξ = (ξ 1 , ξ 2 , ξ 3 ). 2. Multiply the frequency window U j,l (ξ) f (ξ) for each scale j and angle l. 3. Wrap around the origin and obtain W( U j,l f )(ξ), where the range for ξ p is -L p,j,l /2 ≤ ξ p < L 1,j,l /2, j = (0, j f ). No wrapping at scales 0 and j f . 4. Apply 3D inverse FFT to each W( U j,l f ) to obtain the discrete coefficients c D µ .

The computational complexity of the DCuT is O(n 3 log n) flops for n × n × n data [START_REF] Candès | Fast discrete curvelet transforms[END_REF].

INRIA

Reconstruction of threshold-processed curvelet coefficients

u c = µ τ (c D µ (f )) ϕ µ , (31)
leads to denoising or extraction of videos. Here τ (x) denotes a thresholding function. Some classical thresholding functions are listed below [START_REF] Antoniadis | Wavelet methods in statistics: some recent developments an their applications[END_REF]: 1) Hard thresholding

τ (x) = x, |x| ≥ σ, 0, |x| < σ, (32)
2) Soft thresholding

τ (x) =    x -σ, x ≥ σ, 0, |x| < σ, x + σ x ≤ -σ, (33)
3) Firm thresholding

τ (x) := x -σ 2 x , |x| ≥ σ, 0, |x| < σ, (34)
4) SCAD by Antoniadis and Fan [START_REF] Antoniadis | Regularization of wavelets approximations[END_REF] τ

(x) =    sign(x) max(0, |x| -σ), |x| ≤ 2σ, (α-1)x-ασsign(x) α-2 , 2σ < |x| < ασ, x, |x| ≥ ασ. (35)
Generally, we choose α = 3.7. Here σ is a given threshold value.

TV-SYNTHESIS 3D CURVELET SHRINK-AGE

It has been well known that tools from computational harmonic analysis suffer from pseudo-Gibbs phenomena (i.e., oscillation artifacts near the discontinuities), although curvelets have much improved the problem in comparison to traditional wavelets. Reconstructing the coefficient using a rule of total variation (TV) minimisation can reduce the pseudo-Gibbs and element-like artifacts [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF][START_REF] Ma | Curvelets for surface characterization[END_REF][START_REF] Ma | Combined complex ridgelet shrinkage and total variation minimization[END_REF]. Following the previous work, here we combine the TV technique with 3D curvelet transform for the first time. The TV-synthesis curvelet shrinkage can reduce the artifacts and sharpen the blurry edges in extracted fields.

For a function u with |∇u| ∈ L 1 (Ω), the total variation functional is defined by [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] T

V (u) = Ω |∇u(x)| dx. (36)
In order to avoid computational difficulties arising from the non-differentiation of the modulus at zero, the TV functional is often replaced by

T V (u) = Ω |∇u(x)| 2 + β 2 dx, ( 37
) RR n°6683
where β > 0 is a small parameter. In the following description, we mainly restrict our attention to the 3D problem. For the 2D problem, we refer readers to [START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF][START_REF] Ma | Curvelets for surface characterization[END_REF], where the authors have combined TV minimisation with the ridgelet and curvelet transform for image processing. Using a first order finite differences scheme for computing the gradient, the total variation functional for u := (u ι,κ,ν, ) (ι,κ,ν)∈I 3 n is given by

T V (u) = ι,κ,ν (δ 1 u) 2 ι,κ,ν + (δ 2 u) 2 ι,κ,ν + (δ 3 u) 2 ι,κ,ν + β 2 dx, (38)
where (δ 1 u) ι,κ,ν = u ι+1,κ,ν -u ι,κ,ν , (δ 2 u) ι,κ,ν = u ι,κ+1,ν -u ι,κ,ν , and (δ 3 u) ι,κ,ν = u ι,κ,ν+1 -u ι,κ,ν . More precisely, for a given u let

U := {u := (u ι,κ,ν ) (ι,κ,ν)∈I 3 n : c D µ> = c D µ , ∀µ ∈ Λ}. ( 39
)
Then we are looking for the solution of the constrained minimisation problem

min u∈U T V (u). ( 40
)
If the linear subspace V consists of functions on I 3 n given by

V := {υ := (υ ι,κ,ν ) (ι,κ,ν)∈I 3 n : c D µ = 0, ∀µ ∈ Λ}, (41)
the idea of TV-minimisation is to remove the pseudo-Gibbs oscillations by minimizing the functional

F (u) = Ω |u -u 0 | 2 dx + λ T V (u) ( 42
)
for u ∈ {u c + v, v ∈ V}, where u 0 is an original flow, u c is a reconstructed flow after curvelet hard thresholding, and V is a linear subspace of functions consisting of the components removed by thresholding. It should be noted that because of the constraint on the subspace V, Eq. ( 42) is not the usual Rudin-Osher-Fatemi's TV model as in Ref. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], but instead is a variant of the TV problem, which was originally inspired by Durand et al. [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF] for wavelets and extended by Ma et al. [START_REF] Ma | Combined complex ridgelet shrinkage and total variation minimization[END_REF][START_REF] Ma | Curvelets for surface characterization[END_REF][START_REF] Ma | Combined curvelet shrinkage and nonlinear anisotropic diffusion[END_REF]. Using u c as an initial guess, the constrained TV-minimisation can be computed by a projected subgradient descent scheme [START_REF] Durand | Reconstruction of wavelet coefficients using total variation minimization[END_REF]

u l+1 = u l -t l P V (g T V (u l )). ( 43
)
Here g T V (u) denotes the subgradient of T V at u. The step size t l can be taken appropriately to ensure convergence. P V (u) denotes a projection of u on the constrained subspace V. This means that only the coefficients with absolute value smaller than a given threshold σ will be changed by the minimisation process. Let T be the curvelet transform and T -1 be its inverse, then we have

P V (u) = T -1 τ -1 T (u)
where τ -1 denotes the so-called inverse thresholding function,

τ -1 (x) := 0, |x| ≥ σ x, |x| < σ. (44)
A crucial step is to compute the gradient of total variation, i.e., g T V (u) or Figure 7 illustrates the 3D grids for computation of ∇ u T V (u) for 3D problems. The four points denoted by red arrow are points that are related to the derivation of the sum with respect to u ι,κ,ν . The arrow head indicates the direction of the Euler forward difference scheme. The derivation of the total variation at location u ι,κ,ν is given by

∇ u T V (u). INRIA u l,κ,ν u l,κ+1,ν u l,κ+1,ν u l+1,κ,ν u l,κ-1,ν u l,κ,ν-1 u l-1,κ,ν
∇ u T V (u) : = (3u ι,κ,ν -u ι+1,κ,ν -u ι,κ,ν+1 -u ι,κ+1,ν ) • A -1/2 +(u ι,κ,ν -u ι,κ-1,ν ) • B -1/2 +(u ι,κ,ν -u ι,κ,ν-1 ) • C -1/2 +(u ι,κ,ν -u ι-1,κ,ν ) • D -1/2 (45)
for the inner points (ι, κ) ∈ I 3 n and corresponding modification at the boundary

∂I 3 n . Here A = (u ι+1,κ,ν -u ι,κ,ν ) 2 + (u ι,κ,ν+1 -u ι,κ,ν ) 2 + (u ι,κ+1,ν -u ι,κ,ν ) 2 , (46)
B = (u ι+1,κ-1,ν -u ι,κ-1,ν ) 2 +(u ι,κ-1,ν+1 -u ι,κ-1,ν ) 2 +(u ι,κ,ν -u ι,κ-1,ν ) 2 , ( 47
) C = (u ι+1,κ,ν-1 -u ι,κ,ν-1 ) 2 +(u ι,κ+1,ν-1 -u ι,κ,ν-1 ) 2 +(u ι,κ,ν -u ι,κ,ν-1 ) 2 , ( 48
) D = (u ι,κ,ν -u ι-1,κ,ν ) 2 +(u ι-1,κ,ν+1 -u ι-1,κ,ν ) 2 +(u ι-1,κ+1,ν -u ι-1,κ,ν ) 2 . (49)
Essentially, TV minimisation does not set the insignificant coefficients to zero as conventional shrinkage does, but typically removes optimally small values to eliminate the artifacts. This 3D TV constrain method can be also incorporated easily into other transforms such as wavelets, contourlets and Surfacelets [START_REF] Lu | Multidimensional Directional Filter Banks and Surfacelets[END_REF].

Numerical experiments

In this section, we will assess abilities of the 3D curvelet transform combined with TV-minimisation based thresholding for denoising and extracting from image sequences some pertinent structures of the observed geophysical fluids.

Simulated fluid motion video data (denoising and edge extraction)

Fig. 8 (a) shows a 3D vortex video of a simulated geophysical flows with 128 frames. It shows the motion of a vortex in a fluid submitted to the Coriolis force. Image pixels correspond to the concentration of a passive tracer which is advected by the velocity field (computed by the model). Fig. 8 (b) shows the result of using reconstruction of detailed coefficients of 3D wavelet transform [START_REF] Daubechies | Ten lectures on wavelets[END_REF]. Fig. 8 (c) is obtained using the 3D curvelet transform described in section 2.

Both transforms can capture the edges. However, curvelet transform results to be smoother and more continues along the edges than wavelet transform does. It can be seen obviously that oscillating artifacts arise along the edges with wavelet transform: it is due to its poor abilities for representing curve-line singularities.

It is also possible to extract edge structures in noisy environment. Fig. 9 shows 3D evolution of a flow with Gaussian white noise. Fig. 9 (b) and (c) show the denoising results by 3D wavelet transform and curvelet transform, respectively. Fig. 9 (d) and (e) are the edge mapping of (b) and (c) using simple reconstruction of detailed coefficients after hard thresholding. The 3D curvelet shows its good performances for edge-preserving denoising and feature extraction. In these above tests, we only use the simplest hard thresholding for wavelet and curvelet transform. Many advance techniques including contrast enhancement (e.g., [START_REF] Starck | Astronomical image representation by the curvelet transform[END_REF]) can be incorporated into this framework. Now we will show the abilities of the proposed 3D TV-synthesis curvelet shrinkage. Figure 10 (a) shows part of the 3D vortex video with white noise. Fig. 10 (b) is the denoising result by curvelet transform. It should be noticed that analysing such kind of images with weak features using conventional methods is a challenging problem. Indeed the edges are too blurry and almost without high-frequency components. Using our TV-synthesis curvelet transform can sharpen the edges to some extent when we remove the noise: Fig. 10 (c). Fig. 10 (d) shows the removed components in our curvelet denoising. Some edge structures have been kicked out yet. There is a big room for future researches in image processing communities. In order to see clear the difference, Fig. 12 shows the close-up of the original flow, wavelet denoising result, and curvelet denoising result. It can be seen obviously that the results by our proposed method preserve the edges much better than those by wavelet method.

The oscillating artifacts and element-like artifacts that appear in wavelet's results are almost suppressed by our TV-synthesis curvelet method. However, the orthogonal wavelet transform is much faster than the curvelet method. Specially for this case, the elapsed time of forward 3D wavelet transform is 0.54 seconds, and inverse transform is 0.45 seconds. The elapsed time of forward curvelet transform and inverse curvelet transform are 1.55 seconds and 1.65 seconds, while the TV-synthesis curvelet transform has to pay 85.0 seconds. Fig. 13 shows edge extraction using wavelet transform, curvelet transform combined with hard thresholding and TV-synthesis curvelet transform. Fig. 14 displays the respective close-up amplification. The curvelet-based methods show good abilities to extract line-like geometric edges. These extracted spatio-temporal geometric structures will be used in our sequent work on image assimilation of geophysical fluids. This is also one of main motivations of this paper.

Video data containing high frequency edges and textures (denoising)

The proposed method can be widely used in other fields. Indeed, the TVsynthesis curvelet transform is more significant for images with high-frequency edges and textures. We give an example of real video data (the data was also used in [START_REF] Lu | Multidimensional Directional Filter Banks and Surfacelets[END_REF]) to further show the outstanding performance of our method, in comparison to existing wavelet methods. Fig. 15 (a) shows the 3D noisy video consisting in line edges and point features. Fig. 15 (b) is the denoising result using 3D wavelet transform. Analogically, Fig. 15 (c) is the edge-preserving denoising result using our TV-synthesis curvelet transform. Fig. 16 displays the close-up of Fig. 15 (a) (b) and (c), respectively. The proposed method preserves the edges better and also achieve much higher signal-noise-ratio SNR = 36.99 than those by wavelet SNR = 34.92 dB. Finally, we should notice that we do not pay attention to criticise the wavelets here but to point out the respective advantages of using wavelets and curvelets. Indeed, wavelets are better at representing point-like features while curvelets are optimal for curve-like edges and structural features. Orthogonal wavelet transform is very fast and rich achievements on wavelet algorithm in other disciplines can be applied in our field. For a 64 3 data, the wavelet forward and inverse transform takes 0.83 seconds and 0.54 seconds while curvelet transform takes 1.29 seconds and 1.27 seconds; For a 128 3 data the wavelet takes 2.29 and 2.17 seconds while curvelet takes 26.82 and 56.26 seconds for the forward and inverse transform. At the present time, a fast parallel algorithm of 3D curvelet transform can be available at www.curelet.org for large-scale computations. Normally, the performances of curvelets for larger scale images are better than smaller scale images: indeed more decomposed levels can be used in large size problems so that the elements of curvelets have more anisotropic directional selectivity. Combining curvelets with wavelets would be promising for general geophysical fluids and natural image sequences.

Conclusion

In this paper, we investigated applications of 3D curvelet transform to extraction of spatio-temporal geometric structures from geophysical fluid flow image sequences and denoising. Extracting such structures and tracking their evolution is essential for future forecast systems. Indeed, new data assimilation techniques (namely Image Sequences Assimilation) which will take into account such spatio-temporal structured data are under development.

A total variation based 3D curvelet reconstruction is proposed in order to sharpen the edges and to suppress oscillating artifacts. It can be applied widely in other fields on image and video processing.

Figure 1 :

 1 Figure 1: Image sequence over Europe provided by the METEOSAT satellite (visible channel, source Météo France).

Figure 2 :

 2 Figure 2: Images of sea surface temperature of the Black Sea provided by the AVHRR satellite (infra-red channel, source NOAA).

Figure 3 :

 3 Figure 3: The elements of wavelets (left) and curvelets on various scales, directions and translations in the spatial domain (right).

Figure 4 :

 4 Figure 4: Frequency domain window U 0 (a) and its support (b).

Figure 6 :

 6 Figure 6: An element of 3D curvelets at a coarse scale (upper row) and fine scale (lower row) is shown in three cross-sections (left column) and isosurface (middle column). The right column shows their frequency support. It can be clearly seen that the element with high resolution in the space domain has low resolution ion in the frequency domain. But the supports in space and frequency domain are both local.

Figure 7 :

 7 Figure 7: Diagram of 3D coordinate grid for computation of ∇ u T V (u). The four points with red arrow are points related to the derivation of the sum with respect to u ι,κ,ν . The arrow head indicates the Euler forward difference scheme.

4. 2 Fig. 11

 211 Fig.11shows denoising result by wavelets and TV-synthesis curvelet transform of a real -then noisy -EUMETSAT (water vapour channel) video data showing the evolution of a vorticity band. Fig.11 (a) is a frame of original flows with noisy. Fig.11 (b) and (c) are the denoising result by 3D wavelet transform and its removed components, respectively. Fig.11 (d) and (e) are results by 3D TV-synthesis curvelet denoising with 0.005 stepsize and 10 iterations, and its removed components. The SNR (signal-to-noise ratio) of wavelet and curvelet's

Figure 8 :

 8 Figure 8: Extraction of edge structure. (a) original vortex flows. (b) extraction by 3D wavelets. (c) extraction by 3D curvelets.

Figure 9 :

 9 Figure 9: Extraction in noisy environment. (a) Noisy flows. (b) denoising by 3D wavelets. (c) denoising by 3D curvelets. (d) extraction by wavelets from (b). (e) extraction by curvelets from (c).

INRIAFigure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :Figure 15 :Figure 16 :

 10111213141516 Figure 10: (a) Noisy flows. (b) denoising by the 3D curvelets. (c) denoising by the TV-synthesis curvelet transform. (d) removed components by (c).

RR n°6683

Acknowledgements

INRIA

The first author would like to thank financial support from NSFC Grant No. 40704019, TBRF (JC2007030), PetroChina Innovation Fund (060511-1-1), and Invited Professorship at Laboratory Jean-Kuntzman of the University Joseph Fourier and INRIA-Grenoble. This work was partially supported by the French Research Agency through the ADDISA project http://addisa.gforge.inria.fr.