Thibaut Archer Abc

Myriam Berny Ac

Pierre Beauchêne

François Hild

Creep behavior identification of an environmental barrier coating using full-field measurements

Keywords: Environmental barrier coating (EBC), Finite Element Model Updating (FEMU), Full-field measurements, Creep

For the use of thermal and environmental barrier coating (T/EBC) with ceramic matrix composites, it is crucial to master the behavior under (extreme) environments representative of the hot section of engine turbines. An experimental setup to simulate such thermal loading has been developed with various diagnostics enabling for kinematic and thermal field measurements, which are used to drive a finite element model and estimate thermomechanical properties such as creep parameters.

Introduction

Environmental barrier coatings (EBCs) have been developed

for challenging environments such as gas turbine engines to improve fuel efficiency [START_REF] Boyce | Gas turbine engineering handbook[END_REF], [START_REF] Walsh | Gas turbine performance[END_REF]. Substantial gains may still be achieved with increasing operating temperatures in the high-pressure turbine section but require new generations of materials to replace the current, internally cooled, superalloy structures. Among various candidates, ceramic matrix composites (CMCs) are currently considered as an attractive option due to their high temperature resistance and low density [START_REF] Evans | The mechanical behavior of ceramic matrix composites[END_REF]- [START_REF] Morscher | Tensile creep and fatigue of Sylramic-iBN melt-infiltrated SiC matrix composites: Retained properties, damage development, and failure mechanisms[END_REF]. However, the presence of water vapor leads to rapid surface recession [START_REF] Opila | SiC Recession Caused by SiO2 Scale Volatility under Combustion Conditions: II, Thermodynamics and Gaseous-Diffusion Model[END_REF]- [START_REF] Verrilli | Effect of Environment on the Stress-Rupture Behavior of a Carbon-Fiber-Reinforced Silicon Carbide Ceramic Matrix Composite[END_REF],

hence the development of protective coatings not only to withstand higher temperatures but also to ensure environmental stability with a low recession rate [START_REF] Lee | Upper Temperature Limit of Environmental Barrier Coatings Based on Mullite and BSAS[END_REF], [START_REF] Lee | Rare earth silicate environmental barrier coatings for SiC/SiC composites and Si3N4 ceramics[END_REF].

Another concern about EBCs is about both quasi-static and time-dependent mechanical properties, which are crucial to understand failure modes and develop reliable systems.

Thermal cycling at high temperatures induces creep and sintering, thereby resulting in through-thickness cracking during cooling [START_REF] Zhu | Influence of High Cycle Thermal Loads on Thermal Fatigue Behavior of Thick Thermal Barrier Coatings[END_REF], [START_REF] Zhu | Determination of creep behavior of thermal barrier coatings under laser imposed high thermal and stress gradient conditions[END_REF], which may accelerate coating failure. It is therefore essential to know the thermomechanical behavior of EBCs. Among existing techniques, Zhu and Miller [START_REF] Zhu | Sintering and creep behavior of plasma-sprayed zirconia-and hafnia-based thermal barrier coatings[END_REF] have determined sintering and low-stress creep parameters using a dilatometer in isothermal conditions on self-supported ceramic coatings.

The results appeared to be similar to those obtained with laser creep/sintering tests [START_REF] Zhu | Determination of creep behavior of thermal barrier coatings under laser imposed high thermal and stress gradient conditions[END_REF], where the system was uniformly heated on the coating surface, and high temperature mechanical tests [START_REF] Thurn | High-temperature deformation of plasmasprayed ZrO2 thermal barrier coatings[END_REF], [START_REF] Johnsen | Compressive Fatigue of a Plasma Sprayed ZrO2-8wt%Y2O3 and ZrO2-10wt%NiCrAICoY TTBC[END_REF].

In all those cases, in order to identify and validate constitutive laws, comparisons between experiments and modeling are required, hence the development of various identification methods [START_REF] Avril | Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements[END_REF], some of them based on full-field measurements [START_REF] Michel Grédiac | Full-field measurements and identification in solid mechanics[END_REF]. Among these techniques, Finite Element Model Updating (FEMU) is widely used [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF]- [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF]. It is an iterative method based on the minimization of the norm of residuals between numerical and experimental data that can be displacement fields (FEMU-U), load levels (FEMU-F) or combinations of both (FEMU-UF) for instance [START_REF] Michel Grédiac | Full-field measurements and identification in solid mechanics[END_REF]. Last, in addition to using measured displacement fields as raw data, some works [START_REF] Gras | Identification of the Out-of-Plane Shear Modulus of a 3D Woven Composite[END_REF]- [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF] have proposed weighted FEMU based on the covariance matrix of measurement errors, which is considered in this study as well.

In the following, a high temperature test on an environmental barrier coated CMC is first presented. Fullfield measurement techniques at high temperature are then introduced with infrared thermography and Digital Image Correlation (DIC). Last, the calibration of creep parameters is proposed using sequentially FEMU based on displacement and temperature fields.

Materials and experimental set-up

Tested materials

The tested sample is a ca. 3 mm thick CERASEP ® A600 SiC/SiC composite made by SAFRAN Ceramics, whose length and width are respectively 95 and 20 mm. The substrate is then coated with a thin layer of (Si) bond coat, and finally with an EBC made of rare earth silicate. The length of the coated surface is 86 mm. The system is then heated in a furnace to stabilize the coating.

As mentioned before, the nonlinear behavior of the coating is essential to understand and quantify failure modes. In order to describe primary creep, a Norton-Bailey law is commonly used [START_REF] Zhu | Influence of High Cycle Thermal Loads on Thermal Fatigue Behavior of Thick Thermal Barrier Coatings[END_REF], [START_REF] Zhu | Determination of creep behavior of thermal barrier coatings under laser imposed high thermal and stress gradient conditions[END_REF], [START_REF] Thurn | High-temperature deformation of plasmasprayed ZrO2 thermal barrier coatings[END_REF] 𝜀̅ ̇𝑐𝑟 = 𝐴𝑟 0 • 𝑒 -𝐸 𝑎 𝑅𝑇 • 𝜎 ̃𝑛 • 𝑡 𝑚 [START_REF] Boyce | Gas turbine engineering handbook[END_REF] with

𝜀̅ ̇𝑐𝑟 = √ 2 3 𝜺̇𝑐 𝑟 : 𝜺̇𝑐 𝑟 (2)
where 𝜀̅ ̇𝑐𝑟 is the equivalent creep strain rate, 𝜺̇𝑐 𝑟 the creep strain tensor, 𝐴𝑟 0 the pre-exponential factor, 𝐸 𝑎 the creep activation energy, 𝑅 the gas constant, 𝑇 the absolute temperature, 𝜎 ̃ the Von mises equivalent stress, 𝑛 the creep stress exponent, 𝑡 time and 𝑚 the creep time exponent.

In the following, the initial value of the creep stress exponent 𝑛 and time exponent 𝑚 were taken from Ref. [START_REF] Richards | Mechanical properties of air plasma sprayed environmental barrier coating (EBC) systems: preliminary assessments[END_REF] for a similar material (i.e., 𝑛 = 0.9 and 𝑚 = 0.78). Over the investigated temperature range, the creep activation energy was provided by Safran and is equal to 272 ± 27 kJ/mol. The test consists of three heating steps, each applied for 1 h. The maximum temperature is limited to 1350°C in order to avoid sintering of the coating, the latter would induce variations in Young's modulus [START_REF] Siebert | Changes in porosity and Young's Modulus due to sintering of plasma sprayed thermal barrier coatings[END_REF]. During the test, thermal and visible light images were continuously acquired with a period of 2s. Figure 2 shows the laser power, dimensionless front and backside temperature profiles.

Experimental set-up and performed test

Full field measurements at high temperature

In this section, measurement of temperature and displacement fields are first detailed in a generic way, and then applied to the test described in the previous section.

3.1 Temperature field measurement

Calibration

Accurate temperature field measurement is essential not only to control the maximum temperature but also to estimate thermal gradients (i.e. the stress level). In quantitative monochromatic thermography, the knowledge of the emissivity is required to assess the true temperature

𝐿 𝜆 (𝑇) = 𝐿 𝜆 0 (𝑇 𝑅 ) = 𝜀 𝜆 (𝑇) 𝐿 𝜆 0 (𝑇) (3)
where 𝐿 𝜆 is the spectral radiance, 𝑇 the true absolute temperature, 𝐿 𝜆 0 the blackbody spectral radiance, 𝑇 𝑅 the radiance or blackbody temperature, and 𝜀 𝜆 the spectral emissivity. As the device used to measure temperature fields is a middle wavelength infrared camera [3 µm ; 5 µm] with a band-pass filter at 4 µm, the emissivity has been measured by CEMHTI (Orléans, France) [START_REF] Massa | High temperature emissivity, reflectivity, and x-ray absorption of BiFeO3[END_REF], [START_REF] Del Campo | High-Temperature Radiative Properties of an Yttria-Stabilized Hafnia Ceramic: High Temperature Radiative Properties of an YSH Ceramic[END_REF] and is shown in The sensor is sensitive to the radiant flux coming from the thermal scene, which is converted into digital levels (DL) and then translated into temperature with the emissivity of the observed surface. Digital levels are expressed with a calibration law [START_REF] Rodiet | Mesure de Température par Méthodes Multi-Spectrales et Caractérisation Thermique de Matériaux Anisotropes par Transformations Intégrales:«Aspects Théoriques et Expérimentaux[END_REF] as

𝑈 = 𝑠 [𝜑 𝑜 (𝑇) + 𝜑 𝑒 (𝑇) + 𝑘(𝑇 𝑐𝑎𝑚 4 -𝑇 𝑑 4 ) ] 𝐼𝑇 + 𝑈 0 ( 4
)
where 𝑈 is the measured tension, 𝑠 the sensitivity of the sensor, 𝜑 𝑜 (𝑇) the radiant flux coming from the object, 𝜑 𝑒 (𝑇)

the radiant flux coming from the environment, 𝑘(𝑇 𝑐𝑎𝑚

4

-𝑇 𝑑 4 ) the radiant flux coming from the camera housing, 𝐼𝑇 the integration time, and 𝑈 0 a digital offset.

The global radiant flux is expressed in Equation ( 4)

as the sum of three contributions. In this case, given the range of temperature, the radiant flux from the object is prevalent. In addition, because of the band-pass filter, the camera is monochromatic and the radiant flux from the object can be directly related with Planck's law for a blackbody. These assumptions lead to some simplifications of Equation ( 4)

𝑈 = 𝑝 • [ 𝐶 1 𝜆 -5 𝑒 ( 𝐶 2 𝜆𝑇 ) -1 ] • 𝐼𝑇 + 𝑈 0 (5)
with 𝑝 and 𝑈 0 constants to identify, 𝐶 1 and 𝐶 2 parameters of Planck's equation. This calibration expression is then identified with the manufacturer data on a blackbody for several integration times as shown in Figure 4. The good agreement enables for the use of Equation ( 5) on the real material.

To check all these assumptions, residuals between a temperature field rebuilt using Equation ( 6) and the builtin calibration data for an emissivity of 0.

Digital image correlation at high temperatures

DIC was used numerous times as a non-contacting full-field measurement technique [START_REF] Michel Grédiac | Full-field measurements and identification in solid mechanics[END_REF] at room temperature [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF]. It has been proven to be efficient in high temperature tests since the pioneering work of Lyons et al. [START_REF] Lyons | Hightemperature deformation measurements using digital-image correlation[END_REF] where the authors tested samples up to 650°C. Tests were then performed at higher temperatures with DIC. Pan et al. [START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200 C[END_REF], Grant et al.

[30] reached temperatures as high as 1200°C, Leplay et al.

reached 1350°C [START_REF] Leplay | Analysis of Asymmetrical Creep of a Ceramic at 1350°C by Digital Image Correlation[END_REF], Novak and Zok [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF] 1500°C. These works highlighted some challenges:

-Gray level conservation may be violated. To circumvent this issue, filtering thermal radiations via lighting with band-pass filters is one possible experimental solution [START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200 C[END_REF], [START_REF] Grant | High-temperature strain field measurement using digital image correlation[END_REF], [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF]. Brighness and contrast corrections during the DIC procedure is another route, be it for local DIC [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications[END_REF] or global approaches [START_REF]Integrated Digital Image Correlation considering gray level and blur variations: Application to distortion measurements of IR camera[END_REF]- [START_REF] Mendoza | Differentiating 3D textile composites: A novel field of application for Digital Volume Correlation[END_REF].

-Fluctuations due to heat haze effects impact the measurement by creating artifacts of the optical path.

Using forced air flux [START_REF] Lyons | Hightemperature deformation measurements using digital-image correlation[END_REF], [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF], or spatiotemporal DIC [START_REF] Besnard | Analysis of image series through global digital image correlation[END_REF], [START_REF] Berny | On the analysis of heat haze effects with spacetime DIC[END_REF] are two possible solutions to mitigate these effects.

-For stability, high temperature speckles are required.

Various coatings made of nitride and aluminum oxides [START_REF] Lyons | Hightemperature deformation measurements using digital-image correlation[END_REF], or cobalt oxide mixed with inorganic adhesive [START_REF] Pan | High-temperature digital image correlation method for full-field deformation measurement at 1200 C[END_REF] have been used. High temperature powders are another potentially durable solution [START_REF] Dong | A Review of Speckle Pattern Fabrication and Assessment for Digital Image Correlation[END_REF], [START_REF] Archer | Global digital image correlation up to very high temperatures with grey level corrections[END_REF].

Global DIC framework

Global DIC considers the minimization of the gray level residuals over the whole region of interest (ROI). It ensures continuity of the displacement fields, which is assumed to be true in this work. The gray level conservation at any pixel location reads

𝑔(𝒙 + 𝒖(𝒙)) = 𝑓(𝒙) (7)
where 𝑓 is the gray level image in the reference configuration, 𝑔 in the deformed configuration, 𝒙 any pixel in the ROI, 𝒖(𝒙) the sought displacement field. The correlation procedure aims to find 𝑢 that minimize the L2norm of the gray level residual 𝜂 defined as

𝜂(𝒙) = 𝑓(𝒙) -𝑔(𝒙 + 𝒖(𝒙)) (8)
As the minimization of the squared norm of 𝜂 over the whole region of interest Ω is an ill-posed problem, the displacement field is decomposed as

𝒖(𝒙) = ∑ 𝑢 𝑖 𝜽 𝒊 (𝒙) 𝒊 (9)
where 𝜽 𝒊 are chosen shape functions and 𝑢 𝑖 the associated nodal displacements to determine. The Gauss-Newton iterative procedure [START_REF] Mathieu | Estimation of Elastoplastic Parameters via Weighted FEMU and Integrated-DIC[END_REF] that solves linear systems written in terms of displacement corrections {𝛿𝒖} reads

[𝑴 𝑫𝑰𝑪 ]{𝜹𝒖} = {𝒃} (10)
where [𝑴 𝑫𝑰𝑪 ] is the Hessian matrix, {𝜹𝒖} the vector gathering all the amplitude corrections and {𝒃} the residual column vector. In this work, mechanical regularization [START_REF] Tomičevć | Mechanics-aided digital image correlation[END_REF] is also used to enforce mechanical admissibility to the displacement fields. Table 3 summarizes the DIC parameters used hereafter. to avoid large gray level changes due to high temperatures [START_REF] Lyons | Hightemperature deformation measurements using digital-image correlation[END_REF], [START_REF] Novak | High-temperature materials testing with full-field strain measurement: Experimental design and practice[END_REF]. This assumption is confirmed in Figure 9 in which the gray level residuals are plotted between the first image and an image in the middle of the last step. It is observed that the mean residual on the black coated surfaces is close to 0 (Root Mean Square residual is equal to 6 GL) whereas dots appear where there are white speckles. This is due to a darkening of the boron nitride coating, which is visible after the test.

Thermal model

The tested sample is modeled as an assembly of two parts, one representing the SiC/SiC composite and one the EBC (Figure 11). Neither the Si bond coat layer nor the grips are accounted for. In the following analysis, the thermomechanical parameters of the CMC part are assumed to be known and will not be calibrated.

where 𝐴 denotes the laser power, 𝐵 𝑥 and 𝐵 𝑦 the widths of the laser beam (respectively, along the 𝑥-axis and 𝑦-axis), 𝑘 the power of the super Gaussian, 𝑋 𝑐 and 𝑌 𝑐 the center coordinates of the laser beam.

The effect of the grips is modeled as heat losses on the top and bottom areas of the CMC material (depicted by red arrows in Figure 11), with a constant amplitude ϕ 0 .

Thermal exchanges are also considered in the heat transfer model. Given the high temperatures, radiation exchanges are prevalent. The total emissivity measured on the EBC and plotted in Figure 3, is used for the top white surface. For all the other surfaces, since they are coated with black paint, the total emissivity provided by the manufacturer is considered. Last, convection is also accounted for with an exchange coefficient depending on the temperature and whether it is the top or bottom surface [43].

Thermomechanical model

The initial hypothesis made on the coating behavior is thermoelasticity, with temperature-dependent thermomechanical parameters. A qualitative comparison is made between measured displacements (elongation and deflection profiles) and computed displacements (Figure 12)

for the first heating step. For all further analyses, the computed displacements (in metric units) are expressed in DIC units (i.e., pixels). experimentally. These variations are mainly due to differences in composition of the processed air plasma sprayed EBC but also in the deposition parameters. Besides, it may also be due to [START_REF] Boyce | Gas turbine engineering handbook[END_REF] an asymmetric creep behavior between tension and compression of the coating, which has already been observed in ceramics [START_REF] Leplay | Analysis of Asymmetrical Creep of a Ceramic at 1350°C by Digital Image Correlation[END_REF], [START_REF] Dusserre | Applicability of image correlation techniques to characterise asymmetric refractory creep during bending tests[END_REF], [START_REF] Massard | Etude du fluage de réfractaires électrofondus du système alumine-zircone-silice[END_REF],

(2) change of creep parameters for higher stress levels as the maximum loading of TMA tests is 40 MPa whereas the compressive stresses reached by the reported test are greater than 150 MPa. Thus, the creep stress and time exponents and the thermomechanical parameters will be updated in the following section.

To complete the thermomechanical model of the experiment, the mechanical boundary conditions are assessed. In the best case, mechanical boundary conditions would correspond to clamping in the grips.

However, the analysis of displacements close to the grips (left part of the green mesh in Figure 14) shows significant motions, compared to those plotted in Figure 10 in both directions. For this reason, boundary conditions applied to the FE model are those measured, spatially interpolated on the FE mesh and applied to the left vertical line on the red mesh shown in Figure 14. The next section will introduce the identification procedure chosen to assess the parameters of the described thermomechanical model.

Identification procedure

In this section, notations

[𝑿] and {𝑿} refer to spatiotemporal matrix or vector.

Finite Element Model Updating

The FEMU-TU identification procedure [START_REF] Kavanagh | Finite element applications in the characterization of elastic solids[END_REF], [START_REF] Pagnacco | Parameter Identification from Mechanical Field Measurements using Finite Element Model Updating Strategies[END_REF]

𝜒 𝑇 2 ({𝒑}) = 1 𝑁 𝑇 𝑁 𝑡,𝑇 ({𝑻 𝒎 } -{𝑻 𝒄 ({𝒑})}) 𝑡 [𝑪 𝑻 ] -1 ({𝑻 𝒎 } -{𝑻 𝒄 ({𝒑})}) (13)
and

𝜒 𝑈 2 ({𝒑}) = 1 𝑁 𝑈 𝑁 𝑡,𝑈 ({𝑼 𝒎 } -{𝑼 𝒄 ({𝒑})}) 𝑡 [𝑪 𝑼 ] -1 ({𝑼 𝒎 } -{𝑼 𝒄 ({𝒑})}) (14)
The FEMU-TU algorithm iteratively updates the parameters by minimizing thermal and kinematic spatiotemporal residuals (gathered in the weighted right-handside vectors

{𝑩 𝑿 }) {𝜹𝒑} = (𝜔[𝑯 𝑻 ] + (1 -𝜔)[𝑯 𝑼 ]) -1 (𝜔{𝑩 𝑻 } + (1 -𝜔){𝑩 𝑼 }) (15)
The search direction is given by the global Hessian matrix

[𝑯], defined as the weighted sum of thermal and kinematic Hessians

[𝑯] = 𝜔[𝑯 𝑻 ] + (1 -𝜔)[𝑯 𝑼 ] (16)
In a Bayesian approach, the weighting factor between the thermal and kinematic parts (introduced in Equation ( 12))

reads [START_REF] Simoncelli | Bayesian Multi-Scale Differential Optical Flow[END_REF] 𝜔 = 𝑁 𝑇 𝑁 𝑡,𝑇 𝑁 𝑇 𝑁 𝑡,𝑇 + 𝑁 𝑈 𝑁 𝑡,𝑈

A relevant way to assess the sensitivity of the sought parameters is to analyze the Hessian matrices [START_REF] Tarantola | Inverse Problems Theory. Methods for Data Fitting and Model Parameter Estimation[END_REF]. The Hessian matrix of a field 𝑋 depends on the sensitivity matrix [𝑺 𝑿 ] (gathering the sensitivity vectors {𝑺 𝑿,𝒑 𝒊 } of the measured data with respect to an independent variation of each parameter 𝑝 𝑖 )

[𝑯 𝑋 ] = [𝑺 𝑋 ] 𝑡 [𝑪 𝑋 ] -1 [𝑺 𝑋 ] 𝑁 𝑋 𝑁 𝑡,𝑋 (18)
with

{𝑺 𝑿,𝒑 𝒊 } = {𝑿 𝒄 ({𝒑 𝒅𝒊𝒔𝒕,𝒊 })} -{𝑿 𝒄 ({𝒑 𝒓𝒆𝒇 })} 𝜀 (19)
The set of reference parameters is denoted {𝒑 𝒓𝒆𝒇 } and the set of perturbed parameters used to compute {𝑺 𝑿,𝒑 𝒊 } is denoted {𝒑 𝒅𝒊𝒔𝒕,𝒊 } and is equal to the reference set with the 𝑖 th parameter 𝑝 𝑖 varied with an offset 𝜀 (taken here equal to 5% of reference value for all parameters, except for those corresponding to a position where the perturbation is made proportionally to the length of the coating). Thus, the sensitivity analysis consists in computing the solution with the initial set of parameters (reference thermal and kinematic fields) then comparing it with solutions from computations with perturbed parameters.

In the following section, the results of such analysis are presented to assess which parameters of the previously described model (Section 4.1) can be calibrated with regards to uncertainties on thermal and kinematic measurements, and the corresponding identification strategy.

Sensitivity analysis

The following sensitivity analysis is performed by computing only one level of heating close to the first heating step, with both transient and steady states. For the sensitivity analysis, the selected thermal parameters {𝒑 𝒕𝒉 } are the six parameters describing the thermal load ϕ (𝐴, 𝐵 𝑥 , 𝐵 𝑦 , 𝑘, 𝑋 𝑐 , 𝑌 𝑐 ), and the amplitude of bottom surface heat losses ϕ 0 .

Regarding the thermomechanical parameters {𝒑 𝒕𝒉𝒎 }, two material properties are analysed, namely the in-plane thermal expansion 𝛼 𝑥 (equal to 𝛼 𝑦 ) and creep. The first one is parameterized as a first order transformation of the initial law (i.e., 𝛼 𝑥,𝑢𝑝𝑑𝑎𝑡𝑒 = 𝛼 𝑥,𝑖𝑛𝑖𝑡 (𝑎 𝑇 + 𝑏)). The second one is described by the two powers 𝑛 and 𝑚, leading to a total of 12 parameters for the sensitivity analysis.

In the sequel, the effect of the weighting factor 𝜔 between thermal and kinematic data is analyzed. The results are compared to a Bayesian approach, with a value of 𝜔 expressed by Equation ( 17) (and here equal to 0.73, i.e.

the global functional is mostly dependent on thermal data), and for different arbitrary values of 𝜔, prescribed by the user. In the last case, the value of 𝜔 ranges from 0 (pure FEMU-U analysis) to 1 (pure FEMU-T analysis).

The different global Hessian matrices are reported in First, it can be noticed that for all values of the weighting factor 𝜔, the term H 77 is very small, which means that the heat loss parameter ϕ 0 cannot be calibrated.

Conversely, its effect is completely negligible in comparison with all the other parameters. Therefore, it will not be considered in the final identification procedure.

Second, a Bayesian approach for the FEMU-TU functional However, performing a identification through FEMU-U does not appear as a good choice since it will exclude from the analysis 73% of the total amount of data of the experiment (thermal data), leading to higher uncertainties on thermal loading (parameters 𝐴, 𝐵 𝑥 and 𝐵 𝑦 ).

Therefore, it is concluded from this sensitivity analysis that a good identification strategy is to perform two uncoupled FEMU steps. It is beneficial to run first a FEMU-T analysis to calibrate the six thermal parameters (𝐴, 𝐵 𝑥 , 𝐵 𝑦 , 𝑘, 𝑋 𝑐 , 𝑌 𝑐 ), using only IR camera data, then to perform a FEMU-U calibration of the five thermomechanical parameters describing creep and thermal expansion.

Implemented schemes

As concluded for the previous section, the calibration strategy is a two-step uncoupled FEMU-T/FEMU-U procedure.

The algorithm first runs and calibrates the thermal loading parameters {𝒑 𝒕𝒉 } = {𝐴, 𝐵 𝑥 , 𝐵 𝑦 , 𝑘, 𝑋 𝑐 , 𝑌 𝑐 } 𝑡 , for each temperature level, independently, by minimizing the residuals between measured and computed temperatures. Thus, from one step to another, besides its amplitude, the shape of the super-Gaussian (position of laser center, width, and shape) may evolve.

For each step, only the steady state part is considered to get a first (and fast) estimate of the loading parameters.

The algorithm iteratively updates thermal parameters {𝒑 𝒕𝒉 } (according to functional 𝜒 𝑇 2 ({𝒑 𝒕𝒉 } ) by determining the new

correction vector {𝜹𝒑 𝒕𝒉 } {𝜹𝒑 𝒕𝒉 } = [𝑯 𝑻 ] -1 {𝑩 𝑻 } (20)
with [𝑯 𝑻 ] defined in Equation ( 18) and

{𝑩 𝑻 } = 1 𝑁 𝑇 𝑁 𝑡,𝑇 [𝑺 𝑇 ] 𝑡 [𝑪 𝑇 ] -1 ({𝑻 𝒎 } -{𝑻 𝒄 ({𝒑 𝒕𝒉 } )}) (21)
where [𝑪 𝑇 ] is a diagonal matrix with a variance depending on the temperature level to weight temperatures according to their uncertainties. The values of the standard deviation were identified in other tests. The minimization stops once the magnitude of the new correction {𝜹𝒑 𝒕𝒉 } is less than 1%.

Once the thermal loading parameters are calibrated for the three temperature levels, the entire experiment is considered, with transient and steady states, for the FEMU-U identification. Thermal boundary conditions are set and the FEMU-U procedure runs, according to Equation ( 14), to calibrate the thermomechanical parameters {𝒑 𝒕𝒉𝒎 } = {𝛼 𝑥 1 , 𝛼 𝑥 2 , 𝛼 𝑥 3 , 𝑛, 𝑚} 𝑡 by iteratively computing the corrections

{𝜹𝒑 𝒕𝒉𝒎 } {𝜹𝒑 𝒕𝒉𝒎 } = [𝑯 𝑼 ] -1 {𝑩 𝑼 } (22)
with [𝑯 𝑼 ] defined in Equation ( 18) and, considering that

[𝑪 𝑼 ] = 2𝜎 𝑓 2 [𝑴 𝐷𝐼𝐶 ] -1 [48], [49] {𝑩 𝑼 } = 1 2𝜎 𝑓 2 𝑁 𝑈 𝑁 𝑡,𝑈 [𝑺 𝑈 ] 𝑡 [𝑴 𝐷𝐼𝐶 ]({𝑼 𝒎 } -{𝑼 𝒄 ({𝒑 𝒕𝒉𝒎 } )}) (23)
where [𝑴 𝐷𝐼𝐶 ] is a block diagonal matrix constructed with the DIC Hessian matrix [𝑴 𝐷𝐼𝐶 ] introduced in Section 3.2.1, and 𝜎 𝑓 the standard deviation of acquisition noise.

Moreover, to help convergence of each algorithm, an additional Tikhonov regularization is introduced [START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF], with

[𝑯 𝑻𝒊𝒌𝒉,𝑿 ] = 𝜆 𝑇𝑖𝑘ℎ,𝑋 [𝑰] (24)
and

{𝑩 𝑻𝒊𝒌𝒉,𝑿 } = 𝜆 𝑇𝑖𝑘ℎ,𝑋 ({𝒑 𝟎 } -{𝒑})./{𝒑 𝟎 } (25)
where [𝑰] is the identity matrix, {𝒑} the set of current parameters (i.e. {𝒑 𝒕𝒉 } or {𝒑 𝒕𝒉𝒎 } depending on the FEMU-X procedure), {𝒑 𝟎 } the initialization of {𝒑} and 𝜆 𝑇𝑖𝑘ℎ,𝑋 the Tikhonov factor, computed here with respect to the largest eigenvalue 𝜆 𝑚𝑎𝑥,𝑋 of the Hessian matrix [𝑯 𝑋 ] (i.e. 𝜆 𝑇𝑖𝑘ℎ,𝑋 = 𝜆 𝑚𝑎𝑥,𝑋 10 𝑖𝑡𝑖𝑘 ). This regularization leads to a slight modification of each FEMU-X system

{𝜹𝒑} = ([𝑯 𝑿 ] + [𝑯 𝑻𝒊𝒌𝒉,𝑿 ]) -1 ({𝑩 𝑿 } + {𝑩 𝑻𝒊𝒌𝒉,𝑿 }) (26)
The Tikhonov regularization is relaxed step by step, by reducing the value of 𝑖𝑡𝑖𝑘 (i.e., decreasing 𝜆 𝑇𝑖𝑘ℎ,𝑋 ) until it possibly reaches the lower eigenvalues of the Hessian matrix [𝑯 𝑋 ].

The following section presents the results obtained with the proposed identification strategy.

5 Discussion of the results

Thermal identification

The results are first presented and detailed for the first heating level. As mentioned in Section 4.1.1, the parameters to be determined are the laser loading, meaning that only results on the stabilized heating steps are analyzed. Figure 19 shows the Root Mean Square (RMS) residual during the minimization procedure. The RMS residual decreases from 66°C to 20°C, thereby showing the efficiency of the procedure (i.e a three-fold reduction). It is also worth noting that for the first two Tikhonov relaxations, the residuals decrease immediately and are stable. This is not the case for the third one. 𝑚 remains almost constant despite several differences in materials but also testing conditions, namely the loading, which is either tensile or compressive, and the stress level, which is higher using the laser test. Further, the stress exponent 𝑛 calibrated with both FEMU cases tends to increase but remains close to 1, which is usually attributed to diffusion creep in ceramics [START_REF] Chokshi | Characteristics of creep deformation[END_REF], possibly with diffusional mass transport along the grain boundaries [START_REF] Coble | A Model for Boundary Diffusion Controlled Creep in Polycrystalline Materials[END_REF] that could be consistent with the microstructure of the coating. To assess the CTE change, Figure 30 shows the deflection at the end of the sample before and after identification compared to the DIC measurement for the first step. As already shown in Figure 25(a), in which the CTE increased after identification, residuals are slightly higher after the procedure, which was not expected. To explain this result, two hypotheses are proposed. First, the CTE sensitivity is too low to be properly identified with this procedure and this test. As only a small part of the sample is heated and only three heating steps are performed, it is the best case to identify a law written over a very large temperature range, particularly for ceramic materials whose CTE is low. Second, the uncoupled identification procedure may not be optimal, especially for this parameter, which is directly dependent on the temperature field. One possible solution would be to use a FEMU-TU procedure for this parameter and then a FEMU-U for the creep parameters.

Uncertainties on the measured temperature field may also explain this result, especially for lower temperatures for which the IR camera is less sensitive. These measurements were then compared to finite element models of the test to confirm that thermoelasticity was not sufficient to describe the material behavior.

Moreover, the accuracy of the measurements allowed the investigated test to be sensitive to creep as there were large residuals between measured and computed deflections with the initial set of creep parameters.

In order to calibrate these parameters, an uncoupled weighted FEMU-T/FEMU-U framework was followed. The identication strategy was selected using a sensitivity analysis and the study of the Hessian matrix, which revealed to be a valuable tool for test/computation correlations.

Thus, the approach aimed to minimize weighted least squares residuals between measured and computed fields, first regarding temperatures to identify the thermal boundary conditions then kinematic data to assess thermomechanical parameters.

Furthermore, the results after thermal identification showed low residuals on the temperature fields and different results depending on the area of minimization. In both cases, the mechanical identification ended up with creep parameters in which the stress exponent was close to 1, which is consistent with diffusional creep observed in ceramic materials [START_REF] Chokshi | Characteristics of creep deformation[END_REF]. This observation shows the possibility to assess the creep behavior of the coating using full-field measurements combined with laser heating, which has the advantage of being sensitive to a range of temperatures and loadings with only one test. However, because of a globally low sensitivity given the challenging environment, thermomechanical parameters such as the coefficient of thermal expansion were not successfully identified with the chosen procedure and set of parameters.

Concerning the creep parameter identification, other constitutive postulates could be tested. For instance, creep stress and time exponents were considered independent of the temperature. In some cases at high temperatures, this choice may not be valid [START_REF] Thurn | High-temperature deformation of plasmasprayed ZrO2 thermal barrier coatings[END_REF]. Thus, more tests would be required to assess this dependence.

Regarding the full identification procedure, the first step would be to increase the amplitude of displacements.

In that case, it could be performed at the two endpoints of the sample by multiscale image correlation [START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF]. The right side because it is the most sensitive area for both elongation and deflection, and the left side because the kinematic boundary conditions close to the grips should be

Figure 1

 1 Figure 1 displays a high power (3kW) high heat flux CO2 laser locally heating the sample described in Section 2.1. This thermal loading induces multiaxial (in-plane and through-thickness) thermal gradients. Under such conditions, the material expands and deflects. To measure the thermomechanical fields of the sample, the top surface is monitored using a middle-wave infrared (IR) camera (FLIR X6580sc) for temperature fields, a visible light camera (AVT Pike-421 with 4/80 lens and 41 mm extension ring) is positioned to measure the edge motion of the sample thanks to LED lighting (GS VITEC Multiled PT). The backside is monitored using bichromatic (IGAR6, Lumasens) and monochromatic (MI3, Raytek) pyrometers.

FIGURE 1 :

 1 FIGURE 1: High temperature test monitored with infrared and visible light cameras

FIGURE 2 :

 2 FIGURE 2: Laser power in green, dimensionless temperature profiles measured by the infrared camera in red and by the pyrometer in blue during the test

Figure 3 ,

 3 Figure 3, at this wavelength, from a few hundred degrees Celsius to the same maximum temperature as in Figure 2. The measured total emissivity is also plotted and will be used in Section 4.1.1.

FIGURE 3 :

 3 FIGURE 3: Spectral emissivity at 4 µm and total emissivity measured on the environmental barrier coating

FIGURE 4 :

 4 FIGURE 4: Calibration data by FLIR ® manufacturer and estimation for several integration times

  5, and an integration time of 25 µs are shown in Figure 5(b) for a heating with high thermal gradients (Figure 5(a)). heterogeneous residuals are observed with maximum values at both highest and lowest DLs, which is explained for low temperatures by getting close to the lower bound of the calibration range.

FIGURE 5 :Figure 6

 56 FIGURE 5: (a) DL field measured on the studied EBC, and (b) residuals (°C) between Equation (6) and the manufacturer calibration data

FIGURE 6 :

 6 FIGURE 6: Average dimensionless temperature fields at the three stabilized heating steps

Figure 7

 7 Figure7shows the reference image of the monitored edge of the sample during the test. To get a homogeneous background, a high temperature black paint (AREMCO ® 840-M) was applied first, followed by a boron nitride white paint to obtain a random pattern. During the test, the left part of the sample was held in grips.

FIGURE 7 :FIGURE 8 :

 78 FIGURE 7: Reference image of the coated CMC with its speckle pattern. The red dot on the right side of the sample corresponds to the location of the point for which temporal profiles are plotted hereafter.

FIGURE 9 :

 9 FIGURE 9: Gray level residuals (GL) between the first image and an image in the middle of the last step

FIGURE 10 :

 10 FIGURE 10: Displacements of the red dot (Figure 7) during the test (1 pixel is equal to 55 µm)

FIGURE 11 :

 11 FIGURE 11: Thermomechanical model with thermal loads (green arrows) from the laser heating and heat loss in the grips (red arrows)

FIGURE 12 :FIGURE 13 :

 1213 FIGURE 12: (a) Elongation (in blue) and (b) deflection (in orange) at the end of the sample obtained by DIC (solid line) and thermoelastic computation (dashed line), 1 pixel is equal to 55 µm

FIGURE 14 :

 14 FIGURE 14: Focus on the area close to the grips with both DIC (in green) and FE (in red) meshes

  aims to calibrate the sought parameters {𝒑} of an FE model by minimizing the global functional 𝜒 𝑇𝑈 2 on thermal (𝑇) and displacement (𝑈) fields𝜒 𝑇𝑈 2 ({𝒑}) = 𝜔 𝜒 𝑇 2 ({𝒑}) + (1 -𝜔)𝜒 𝑈 2 ({𝒑})(12)where 𝜒 𝑇 is the thermal cost function, 𝜒 𝑈 the kinematic cost function, and 𝜔 the weighting factor between thermal and kinematic contributions.The FEMU-TU procedure relies on a simultaneously weighted[START_REF] Gras | Identification of the Out-of-Plane Shear Modulus of a 3D Woven Composite[END_REF]-[START_REF] Gras | Identification of a set of macroscopic elastic parameters in a 3D woven composite: Uncertainty analysis and regularization[END_REF] least squares minimization of thermal and displacement residuals. The residuals for a field 𝑋 (where 𝑋 denotes the temperatures 𝑇 or nodal displacements 𝑼) is defined as the difference between the measured data {𝑿 𝒎 } and the computed ones {𝑿 𝒄 ({𝒑})}, where the latter ones are projected in space on the measurement mesh and in time on the experimental time basis. The weighting of each cost function is made by considering the number of data acquired and used for the minimization (in space 𝑁 𝑋 and in time 𝑁 𝑡,𝑋 ), and by the associated uncertainties on measurements of the field 𝑋 (through the covariance matrix [𝑪 𝑿 ] discussed in the following section). Therefore, the thermal and kinematic cost functions are written as

Figure 15 .

 15 Figure 15. The Hessian is represented as a square matrix of size the number of parameters, where the values of the absolute components are expressed in logarithmic scale (base 10) and where the sign of the components (plus or minus) of the matrix is plotted for each component. Thus, a high value of the component of the diagonal of the Hessian matrix corresponds to a high sensitivity of such parameter.Couplings between parameters are quantified by the offdiagonal components and can also be assessed with the eigenvectors of the diagonalized Hessian matrix.

Figure 15 :

 15 Figure 15: Global Hessian matrix in log-10 scale with different weighting factors 𝜔

  (i.e 𝜔 = 0.73) leads to a high sensitivity of the thermal parameters and makes the identification of the thermomechanical properties of the coating more difficult. The analysis of the eigenvalues and eigenvectors of the Bayesian Hessian matrix (Figure 16) highlights its bad conditioning. It can be seen that there are 3 to 4 orders of magnitude between the first eigenvalues (whose associated eigenvectors are driven by thermal loading parameters) and the eigenvalues associated with creep and thermal expansion contributions.

FIGURE 16 :

 16 FIGURE 16: Dimensionless (a) eigenvalues and (b) eigenvectors of the global Bayesian Hessian matrix (𝜔 = 0.73)

FIGURE 17 :

 17 FIGURE 17: Dimensionless (a) eigenvalues and (b) eigenvectors of the global matrix computed with 𝜔 = 0.25

FIGURE 18 :

 18 FIGURE 18: Dimensionless (a) eigenvalues and (b) eigenvectors of global Hessian matrix computed with 𝜔 = 0

FIGURE 19 :

 19 FIGURE 19: Root mean square residual (°C) history over the entire sample for the first heating step

FIGURE 20 :

 20 FIGURE 20: Change of (a) Y position and (b) the amplitude of the laser loading during the iterative procedure for the first loading step

FIGURE 21 :

 21 FIGURE 21: Residual fields (°C) considering minimization on the entire sample for (a) the first, (b) second and (c) last step. The blue box shows the central area where a second FEMU analysis will be run

FIGURE 22 :FIGURE 23 :

 2223 FIGURE 22: Residual fields (°C) considering the minimization on the central area (see Figure 21(b)) for (a) the second and (b) last steps

FIGURE 24 :

 24 FIGURE 24: Root mean square (a) displacement residual (px) and (b) increment of parameter history during the minimization procedure

FIGURE 25 :

 25 FIGURE 25: Change of (a) coefficient of thermal expansion law and (b) creep stress exponent 𝒏 with the iteration number and Tikhonov relaxation

FIGURE 26 :

 26 FIGURE 26: Root mean square (a) displacement residual (px) and (b) increment of parameter history

FIGURE 27 :

 27 FIGURE 27: Change of (a) coefficient of thermal expansion law and (b) creep law stress exponent 𝒏

Figure 28

 28 Figure 28 shows the deflection at the end of the sample before and after identification compared to DIC measurements. The updated creep parameters reduce the gap between measurements and modeling, with different values as those calibrated during TMA tests. The small remaining difference concerns the beginning of the step, when the temperature field is not stabilized. It can be linked with temperature differences in the transient stage. However, the thermal calibration procedure could also take into account material properties such as diffusivity, which are relevant when trying to reduce the gap in transient stages.

FIGURE 28 :

 28 FIGURE 28: Deflection at the end of the sample obtained by DIC (solid line), for the initial and calibrated computation

FIGURE 29 :

 29 FIGURE 29: Residual displacement fields after identification in the transverse direction for three stabilized heating steps

FIGURE 30 :

 30 FIGURE 30: Elongation at the end of the sample obtained by DIC (solid line), for the initial and updated computation

Table

TABLE 1 :

 1 DIC hardware parameters

		AVT PIKE F-421 (monochrome
	Camera
		detector)
	Camera definition	2048 x 2048 pixels
		Schneider Kreuznach Band pass
	Filter
		Filter, BP 540-70 HT
	Grey Levels
		8 bits
	amplitude
	Lens	Schneider Kreuznach 2.8/50
	Field of view	see text
	Image scale	see text
	Image acquisition
		0.5 fps
	rate
	Exposure time	25 ms

TABLE 2 :

 2 IR hardware parameters

	Camera	FLIR X6580SC
	Camera definition	640 x 512 pixels
	Digital Levels
		14 bits
	amplitude
	Field of view	see text
	Pixel size	200 µm
	Image acquisition
		0.5 fps
	rate
	Integration time	30 µs

TABLE 3 :

 3 DIC analysed parameters

	DIC software	Correli 3.0 [42]
	Image filtering	none
	Element length	20 px
	Regularization length	200 px
	Shape functions	linear (T3)
	Mesh	regular
	Interpolant	cubic
	Displacement noise-
		0.02 px
	floor
	3.2.2 Application to the test

TABLE 4 :

 4 Creep parameters before and after the FEMU-T, FEMU-U procedures

			FEMU-T on	FEMU-T on
	Creep
		Initialization	the entire	the central
	parameters
			surface	zone
	n	0.9	1.32	1.22
	m	0.78	0.8	0.79

measured as accurately as possible. Another possible route to improve the calibration could be to perform more tests with only one heating step at the same temperature but with several laser beam shapes to assess the effect of thermal loading and temperature level separately on the creep properties. The sensitivity analysis performed herein would help to design the best (i.e., most sensitive) loading history.

Acknowledgement

The authors would like to thank Dr. Manon Fernandez for her help and advice in this work, which was supported under PRC MECACOMP, a French research project co-funded by DGAC and SAFRAN Group, managed by SAFRAN Group and involving SAFRAN Group, ONERA and CNRS.


M A Arshad

J Asseline

D Blavet

G De Noni

J Y Laurent

J C Leprun

SOIL PROPERTIES AS AFFECTED BY DIFFERENT LAND USE PRACTICES IN THE LANGUEDOC REGION OF SOUTHERN FRANCE

INTRODUCTION :

The Languedoc-Roussillon region of southern France, mostly dominated by vineyard monoculture, consists of three principal landscapes: scrubland (garrigue), fallow (jachère) and vineyards. This region accounts for approximately 30% of the total area of about 954 740 ha in vineyard production in France. The vineyards (330 335 ha, 1992 data) also accounts for 29.5 % of the area under fruit crops. The intense rain storms, sloping relief, low soil fertility and the ancient human occupation of the mediterranean area along with the recent mechanized farming systems make the soil susceptable to erosion. Intensive mechanical cultivation of soil under vineyard production results in soil erosion. Heavy runoff transports fine soil particles, large quantities of nutrients and pesticides into water streams adversely affecting the quality of soil and water. For the past several years, the vineyard growers have been worried by the continuing fall in production, considered from 20 to 30% in 1998 (Le Monde, 12 September, 1998). The decline in production is attributed to the land degradation resulting from the modern mechanized forms of cultivation. However, the data on the status of the soil quality is lacking. This study was therefore initiated in an area representative of the Languedoc-Rousillon region. Our objective was to determine changes in surface soil properties resulting from three different land use practices in the Languedoc region of southern France.

MATERIAL AND METHODS :

The study site is located about 1 km south of Corconnes town (~ 30 km northwest of Montpellier). The area is characterised by 3 principal types of land use : the vineyard, scrubland and fallow . The climate is Mediterranean, with an annual precipitation of 950 mm (1951-1980 average). The most rainy months are September, October and November, and the rainfall intensity varies from 55mm/hr (once in 10 years) to 80mm/hr (once in 100 years). The soils in the area are classified as Calcic Luvisol (Haploxeralf -US System) and are silt loam in texture, developed on limestone/ blue marls jurassic formation. While the soil under garrigue and fallow have well developed A and B horizons, the A horizons of the vineyard soils are completely removed. The land use systems compared with each other were : undisturbed scrubland (garrigue -G), predominantly oak, thyme, rosemary and grasses (Brachypodium retusum), 25-yr old fallow (F) with dense grass cover (Brachypodium retusum), young vineyard for 5-yr (YV) with inter-row spaces bare of vegetation, controlled by glyphsate and vineyard for 15-yr (AV) mulched with dry vine shoots. It is a common practice that the growers prepare their fields by using heavy machinery (e.g. bulldozers) scraping the A horizons. This is also true for the two vineyards soils selected for this study. Five replicates from each treatment were sampled to determine soil properties ; rainfall simulation and related measurements were done in triplicate. Soil samples were taken from 0-5 cm depth. Bulk density was estimated by determining the oven dry (105 °C) weight of a soil core of known volume. Carbon (after removal of carbonates by acidification) and nitrogen were determined by NA2000 Analyzer, pH by Zeromatic pH metre, and texture and aggregate-size distribution by wet sieving and laser diffraction technique (Malvern instruments, UK). The rainfall simulation study was conducted in triplicate on March 2-5, 1999 with the ORSTOM type apparatus [START_REF] Asseline | Construction et mise au point d'un infiltromètre à aspersion[END_REF]. The tests were carried out on an area of 1m 2 , enclosed by a metal frame sunk into the soil. Each rainfall simulation was conducted for one hour with an intensity of 60mm/hr, representing an event with a return period of >10 yr. Runoff was collected and measured at 1-min. intervals to the total time of 1 hr in a series of bottles at the downslope end of the frames. The antecedent (initial) and saturation(final) soil moisture at 0-5cm depth was determined gravimeterically. Statistical analyses were done by using Statistica software and Newman and Keuls tests.

RESULTS AND DISCUSSION :

The soils under the three ecosystems are silt loam in texture (~ 50% silt, 20% clay and 30% sand). The A1 horizons of soils under garrigue (G) and fallow (F) have well-developed friable and granular structure with low bulk density (~ 1.16 g/cc for G and 1.19 g/cc for J), and high porosity and biological (roots and fauna) activity. The G soil also has a 5 cm thick compact layer just below the A1 horizon. The vineyard soils are more compact (bulk density ~ 1.47 g/cc) which is attributed to the loss of the organic-rich mineral horizons and the influence of the heavy machinery used in initial clearing of the land and the subsequent vineyard production for several years. Despite the textural similarity, soils under the G and F had 6 to 7 times more macroaggregates (1-2 mm) than soils under the vineyards. All the ecosystems differed significantly from each other in organic C and N. Undisturbed soil under G had the maximum C (41.8 mg/g) followed by the F system (30.5 mg/g) while the intensely cultivated soil (YV) had the minimum C contents (5.4 mg/g). The loss of C, presumably by mineralization, was the highest in the YV soil, as indicated by its low C:N (11.7 in YV vs~ 17.2 in G). The data showed much higher amounts of C under the natural garrigue and grassland fallow systems than in vineyard soils. We attribute this high level of C to the addition of freshly decomposed residue [START_REF] Angers | Recently deposited organic matter in soil water stable aggregates[END_REF][START_REF] Puget | Total and young organic matter distribution in aggregates of silty cultivated soils[END_REF]. In our rainfall simulation study, the maximum runoff was observed in the garrigue and vineyard ecosystems, although the highest soil losses occurred in the vineyard soils. No significant soil loss occurred in the G and F soils because of better protective cover and soil structure. The soil loss for the G, F, YV and AV ecosystems were 3.89, 0.00, 160.3 and 83.3 g/m2, respectively. Lower soil loss in the AV than in the YV was attributed to the protective cover created by the mulching practice adopted by the growers. Incorporation of mulch in the AV system also increased the organic matter in this soil, compared to the YV. No runoff was observed in the F while 50% of the rainfall was recorded as runoff for the G system. The two vineyard soils retained the least amount of rain ('imbibition' water; 9 mm) followed by the G soil (23 mm) and the F system (42 mm). In summary, our data demonstrate that the inappropriate management systems currently in use for vineyard production at the study site resulted in destruction of soil structure, loss of organic matter and deterioration of environmental quality.


Michel Narozny

Michel Barret

Dinh-Tuan Pham
email: dinh-tuan.pham@imag.fr

ICA based algorithms for computing optimal 1-D linear block transforms in variable high-rate source coding

teaching and research institutions in France or abroad, or from public or private research centers.

ICA based algorithms for computing optimal 1-D linear block transforms in variable high-rate source coding

Introduction

Transform coding has been extensively developed for coding natural signals like sounds, images (e.g., JPEG [START_REF] Wallace | Overview of JPEG (ISO/CCITT) still image compression standard[END_REF]) and video. The motivating principle of transform coding is that simple coding may be more effective in the transform domain than in the original signal space. Generally "simple coding" corresponds to the use of scalar quantization and scalar entropy coding, because they provide a good trade-off between computational complexity and performance [START_REF] Goyal | Transform coding with backward adaptive updates[END_REF].

In a conventional transform coder [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF][START_REF] Goyal | Theoretical foudations of transform coding[END_REF], an input vector X is first transformed into another vector Y = TX of the same dimension. The components of that vector are then described to the decoder using independent scalar quantizers on the coefficients. Finally, the decoder reconstructs the quantized transform vector Y and then uses a linear transformation U to get an estimate of the original input vector X = U Y. In this paper, we consider mean-square distortion : E[ X -X 2 ], where E denotes mathematical expectation. When the input vector is Gaussian, the optimal transforms T and U satisfy both the conditions T is an orthogonal matrix (i.e., T -1 = T T , where T denotes transposition) that produces uncorrelated transform coefficients and U = T -1 . Such a transform T is often called Karhunen-Loève (KLT) in coding lessons.

The optimality of the KLT and its inverse is well known for high rates [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF], or when optimal fixedrate quantizers are employed [START_REF] Huang | Block quantization of correlated Gaussian random variables[END_REF], or more generally when the quantizers are scale-invariant [START_REF] Goyal | Transform coding with backward adaptive updates[END_REF] for both the fixed-rate and the variable-rate coding models. Moreover, for non Gaussian sources, it is well known that the KLT can be suboptimal in transform coding [START_REF] Effros | Suboptimality of the Karhunen-Loève transform for transform coding[END_REF]. Research in transform coding theory, including subband coding, has been constant for several decades. There is a great interest in extending theoretical results on transform optimality towards low bit-rates [START_REF] Mallat | Analysis of Low Bit Rate Image Transform Coding[END_REF] [START_REF] Goyal | Transform coding with backward adaptive updates[END_REF], or in extending the KLT optimality in more general situations by weakening the assumptions that ensure the optimality. In [START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF],

it is shown that under the classical assumptions of high-rate quantization, mean-square distortion and variable-rate coding, the KLT remains optimal for a more general class of signals than Gaussian ones.

In [START_REF] Mary | A theoretical high-rate analysis of causal versus unitary online transform coding[END_REF], the quantization effect on backward adaptive transform coding of Gaussian sources is studied; the perturbation due to quantization on coding gain is computed up to second order for KLT and causal unitary transform.

Furthermore, under the high resolution quantization hypothesis, nearly everything is known about the performance of a transform coding with entropy constrained scalar quantization and mean square distortion. It is then straightforward to find a criterion that, when minimized, gives the optimal linear transform under the above mentioned conditions. Nevertheless, the optimal transform computation is generally considered as a difficult task [START_REF] Mallat | Analysis of Low Bit Rate Image Transform Coding[END_REF] and the Gaussian assumption is then used in order to simplify the calculus. In this paper we resolve the problem of computing the optimal transform under high-rate entropy constraint scalar quantization hypothesis.

In 1963, Huang and Schultheiss [START_REF] Huang | Block quantization of correlated Gaussian random variables[END_REF] showed that, for scalar Lloyd-Max quantizers, if the vector source X is Gaussian then the mean square distortion is minimized by choosing U = T -1 when T is a variable Y , a Lloyd-Max scalar quantizer with cells S 1 , . . . , S k and outputs y 1 , . . . , y k satisfies the centroid condition y i = E[Y | Y ∈ S i ] for (1 ≤ i ≤ k); and in the other hand the following property satisfied by Gaussian sources when T is a KLT: the jth component of the quantization noise Y -Y and the ith component of Y are uncorrelated for i = j. When the Gaussian assumption does not hold, it is clear that the previous property is satisfied if the components of Y are statistically independent (indeed Y j -Y j and Y i are deterministic functions of Y j and Y i respectively). Hence their proof remains valid under this condition. Moreover, the result by Huang and Schulteiss still holds when the Lloyd-Max scalar quantizers are replaced by uniform scalar quantizers, under the additional assumption of high-rate quantization [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] (since they satisfy the centroid condition). In this paper, we suppose the transform used by the decoder satisfies U = T -1 . This assumption is partially justified by the fact that the linear transforms T we consider in the following give transformed vectors Y with minimal mutual information between components, and hence with components generally close to independence. However, the total independence is rarely achieved and it will be interesting to complete our study by the general case of arbitrary invertible matrices T and U.

In this paper, we first show that under the high-rate entropy constraint scalar quantization hypothesis, for mean-square distortion and the condition U = T -1 , the criterion that gives -when minimized -the optimal linear transform T can be expressed as a classical contrast of Independent Component Analysis (ICA) (actually the opposite of such a contrast), modified by an additional term which can be explained as a pseudo distance to orthogonality. This new presentation of a classical result gives an interesting point of view covering ICA and data compression. Indeed, it is well known that images of natural scenes are well modeled with ICA [START_REF] Bell | The 'independent components' of natural scenes are edge filters[END_REF], this underlies the potential usefulness of ICA to compression [START_REF] Puga | An experiment on comparing PCA and ICA in classical transform image coding[END_REF][START_REF] Marusic | ICA-FIR based image redundancy reduction[END_REF][START_REF] Ferreira | Class-adapted image compression using independent component analysis[END_REF]. Then, we present two variants of the icainf algorithm by Pham [14] that permit to compute the optimal transform T with 1) no constraint and 2) the orthogonality constraint.

Finally, experimental results on both synthetic data and natural images are given in Section 4, they show that the transforms returned by the new algorithms can 1) outperform the KLT when used with synthetic signals, and 2) achieve slightly better performance for high-rates and better visual quality (preservation of lines and contours) than the KLT and 2D DCT when applied to the medium-to-low bit rate compression of natural images. This last result is unexpected, since no optimality is ensured at low bit-rates. The paper begins with a brief review of high bit-rate transform coding. The results presented here have been partially published in [START_REF] Narozny | Modified ICA algorithms for finding optimal transforms in transform coding[END_REF] and [START_REF] Narozny | ICA-based algorithms applied to image coding[END_REF].

High bit-rate transform coding

The general structure of a transform coding scheme is shown in Fig. 1. The class of signals to be encoded is represented by a random vector X = [X 1 , . . . , X N ] T of size N . Although these signals may be multidimensional like images, they are indexed by an integer to simplify notations: X(m). The components of X(m) are successive samples from a source signal (x(n)) n . A conventional transform coder applies a linear invertible transform T : R N → R N to X in order to obtain a random vector Y = [Y 1 , . . . , Y N ] T better suited to coding than X. To construct a finite code, each coefficient Y i is first approximated by a quantized value Y i . We concentrate only on scalar quantizers. The quantized coefficients are then entropy coded. The coded representation is stored or communicated over an errorcorrected (lossless) channel. The receiver (decoder) provides an approximation X = [ X 1 , . . . , X N ] T of the original signal X by applying a linear transform U to the quantized signal Y. In this paper we assume U = T -1 .

A relevant problem for a system designer is to minimize the reconstruction error under bit rate constraint. Here, the optimization criterion is the mean-square distortion

D = 1 N E[ X -X 2 ] which satisfies D = 1 N N i=1 w i E (Y i -Y i ) 2 + 1 N N i=2 i-1 j=1 w ij E (Y i -Y i )(Y j -Y j ) , (1)
where [A] i,j denotes the element localized on the i-th row and the j-th column of A,

w i = N j=1 [T -1 ] 2 j,i (1 ≤ i ≤ N ) and w ij = N k=1 2[T -1 ] k,i [T -1 ] k,j (1 ≤ j < i ≤ N ).
The last term in (1) vanishes when 1) w ij = 0 for i = j, i.e., the T transform is orthogonal -or more generally when the column vectors of T -1 are pairwise orthogonal -or 2) when the quantization noises of different quantizers are uncorrelated and centered. The last condition 2) is satisfied when i) the vector X is Gaussian and the transform T is a KLT 1 or when ii) the transform coefficients are statistically independent -the Gaussian assumption is then useless -and the quantizers satisfy the centroid condition. Moreover, experimental tests show that condition 2) approximatively holds under high-rate quantization hypothesis for any kind of signal.

In the following, we assume that the end-to-end distortion D can be well approximated by a weighting sum of the distortion

D i = E[(Y i -Y i ) 2 ] of each transform coefficient as D ≈ 1 N N i=1 w i D i . (2)
Approximation ( 2) is an equality when T is orthogonal or when Y has independent components; it is a good approximation for any distortion when the transform coefficients are close to independence and last it is a good approximation for any T under the high-rate quantization hypothesis.

Entropy-constrained scalar quantization

A scalar quantizer Q approximates a random variable Y by a quantized variable

Y = Q(Y ). Q is a mapping from a source alphabet R to a reproduction codebook C = { y i } i∈K ⊂ R,
where K is an arbitrary countable index set. We denote p i = P { Y = y i }. The Shannon theorem [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] proves that the entropy H( Y ) =i p i log 2 p i is a lower bound of the average number of bits per symbol used to encode the values of Y . Arithmetic entropy coding [START_REF] Cover | Elements of Information Theory[END_REF][START_REF] Witten | Arithmetic coding for data compression[END_REF] achieves an average bit rate that can be arbitrarily close to the entropy lower bound; therefore, we shall consider that this lower bound

D ≈ cσ 2 2 -2R , with c = 2 2h( Ỹ ) 12 , (3)
where the constant c depends only on the pdf shape.

Optimal bit allocation

Coding (quantizing and entropy coding) each transform coefficient Y i separately splits the total number of bits among the transform coefficients in some manner. This bit allocation problem can be stated this way: one is given a set of quantizers described by their distortion-rate performances as

D i ≈ c i σ 2 i 2 -2R i , R i ∈ R i for (1 ≤ i ≤ N ).
Each set of available rates R i is a subset of the nonnegative real numbers and may be discrete or continuous. The problem is to minimize the end-to-end distortion D in eq. ( 2)

given a maximum average rate R = N -1 N i=1 R i .
It results of the mean theorem (i.e., the arithmetic mean of the w i D i s is not smaller than their geometric mean, with equality if and only if all the terms are equal) that, under the constraint of a given average rate R, the distortion D is minimum if and only if all the w i D i s are equal, in which case

D T (R) ≈ N i=1 w i c i 1/N N i=1 σ 2 i 1/N 2 -2R , (4)
where σ 2 i is the variance of Y i and c i is the constant associated with the standardized variable of Y i according to the relation [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. When no transform T is applied to X, or equivalently when T is the identity I, the minimum distortion associated with the same maximum average rate R is given by

D I (R) ≈ N i=1 c i 1/N N i=1 σ 2 i 1/N 2 -2R
, where σ 2 i is the variance of X i , and c i is the constant associated with the standardized variable of X i according to the relation (3).

Generalized coding gain and maximum reducible bits

In this paragraph we present a criterion called the generalized coding gain (resp. the generalized maximum reducible bits) which is a generalization of the coding gain [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] (resp. the maximum reducible bits) to non-Gaussian signals and non orthogonal linear transforms.

The distortion rate (4) can be used to define a figure of merit that we call the generalized coding gain

G = D I (R) D T (R) ≈ N i=1 c i 1/N N i=1 σ 2 i 1/N N i=1 w i c i 1/N N i=1 σ 2 i 1/N . (5)
It is the factor by which the distortion is reduced because of the linear transform T, assuming high rate quantization and optimal bit allocation. Taking the inverse of (4), we obtain the optimal rate distortion: R T

(D) ≈ 1 N N i=1 h(Y i ) + 1 2N N i=1 log 2 w i -1 2 log 2 [12D]
, from which we can define the

generalized maximum reducible bits R = R I (D) -R T (D) = 1 2 log 2 G
. It is the quantity by which the average number of bits to code X is reduced because of the transform T.

If we assume that X is Gaussian, then Y is Gaussian and we have

c i = c i = πe 6 (1 ≤ i ≤ N ), hence the generalized coding gain becomes G = N i=1 σ 2 i 1/N N i=1 w i σ 2 i 1/N . Furthermore, if
we suppose that T is orthogonal, then the w i s are all equal to one and the generalized coding gain is identical to the well-known coding gain G appearing in texts and scholarly books, which is maximized for any Karhunen-Loève basis of X (see, for example, [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]).

8

Using the following relations (see, e.g., [START_REF] Cover | Elements of Information Theory[END_REF] for notions of information theory)

h(X) = N i=1 h(X i ) - I(X 1 ; . . . ; X N ), h(Y) = N i=1 h(Y i ) -I(Y 1 ; . . . ; Y N ), h(Y) = h(X) + log 2 | det T|, where I(X 1 ; . . . ; X N )
and I(Y 1 ; . . . ; Y N ) denote respectively the mutual information between the components of X and of Y, h(X) and h(Y) the differential entropy of vectors X and Y respectively, the generalized maximum reducible bits can be expressed as follows:

R = 1 N I(X 1 ; . . . ; X N ) - 1 N I(Y 1 ; . . . ; Y N ) - 1 N log 2 | det T| - 1 2N log 2 N i=1 w i . (6)
Moreover, according to the expression of the w i s, and noting Diag(C) for the diagonal matrix having the same main diagonal as C, the last two terms in ( 6) are equal to

-1 2N log 2 det[Diag(T -T T -1 )] det[T -T T -1 ]
.

Hence we deduce the following expression of the generalized maximum reducible bits:

R = 1 N I(X 1 ; . . . ; X N ) - 1 N I(Y 1 ; . . . ; Y N ) - 1 2N log 2 det[Diag(T -T T -1 )] det[T -T T -1 ] . (7)
Let us remark, according to Hadamard's inequality (see e.g., [START_REF] Cover | Elements of Information Theory[END_REF]), that for any definite positive matrix

C of order N , the quantitity log 2 det[Diag(C)] det C
is greater or equal than zero, with equality if and only if the matrix C is diagonal.

It is now clear that the problem of finding a linear transform T which maximizes the generalized coding gain G defined in ( 5) is the same problem as finding T which maximizes R , or equivalently, finding the linear transform T which minimizes the contrast

C(T) = I(Y 1 ; . . . ; Y N ) + 1 2 log 2 det[Diag(T -T T -1 )] det[T -T T -1 ] . (8)
The first term of ( 8) is a measure of the statistical dependence between the transform coefficients Y i .

It is always non-negative, and zero if and only if the variables are statistically independent. As for the second term, it is always non-negative, and zero if and only if the column vectors of T -1 are pairwise orthogonal. Furthermore, if D is a diagonal matrix, one can verify that C(DT) = C(T), i.e., the contrast is scale invariant. Intuitively, this is not surprising since multiplying each component by a factor λ i and each quantization step by the same factor has no impact on both the final rate and the end-to-end distortion. We can now state the following theorem.

Theorem 1 Under the high-rate quantization hypothesis, for entropy-constrained scalar quantizers, mean-square distortion and the condition U = T -1 , a linear transform T is optimal in coding if and only if it minimizes the contrast C(T) of relation [START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF].

Remark now that the contrast C(T) is always non-negative, and that it is equal to zero if and only if

T -1 is a transform with orthogonal columns which produces independent coefficients.

3 Link between ICA and high-rate transform coding

C O (T) = 1 2 log 2 det[Diag(T -T T -1 )] det[T -T T -1 ] . (9)
The second term C O (T) measures a pseudo-distance to orthogonality of the transform T. In general, the optimal transform T opt in transform coding, i.e., the transform which minimizes the contrast (8), will be different from that T ICA which minimizes the first term of (8), i.e., the solution of the ICA problem. It is important to notice here that the classical assumption made in blind source separation problems, that is the observations are obtained from a linear mixing of independent sources, is not really required in the problem of finding the transform that maximizes the generalized coding gain.

The expression of the contrast (8) depends on the definition of the distortion. In this work, we measure the distortion as mean squared error, therefore it is not surprising that orthogonal transforms are favored over other linear transforms since they are energy-preserving.

Modified ICA algorithms for coding

In appendix, we propose two variants of the icainf algorithm by Pham [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF] for the minimization of the contrast [START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF]. The first one, called Generalized Coding Gain Supremum (GCGsup) gives the transform T opt that minimizes the contrast (8), the second, called Orthogonal Independent Composent Analysis (OrthICA) find the orthogonal matrix T orth that minimizes the contrast C ICA (T).

The minimization of the criterion (8) can be done through a gradient descent algorithm, but a much faster method is the Newton algorithm (which amounts to using the natural gradient [START_REF] Amari | Natural gradient works eficiently in learning[END_REF]). As in [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF],

because of the multiplicative structure of our optimization problem, we use multiplicative increment of the parameter T rather than additive increment. Starting with a current estimator T, it consists of expanding C( T + E T) with respect to the matrix E up to second order and then minimizing the resulting quadratic form in E to obtain a new estimate. Note that the parameter E is a matrix of order N . This method requires the computation of the Hessian 2 of C( T + E T) with respect to E, which is quite involved. For this reason, we will approximate it by the Hessian of C( T + E T), computed under the assumption that the transform coefficients Y i are independent. The method is then referred to as quasi-Newton. Although those simplifications result in a slower convergence speed towards the solution, they cause the robustness of the algorithm to be improved by reducing the risk of divergence when the initial estimator T 0 is far from the final solution. Note that the final solution is the same as that obtained without simplification since the algorithm consists of cancelling the first order terms in the expansion of C(T + ET). 2 The Hessian of a function of several variables is the matrix of its second partial derivatives.

11

In this section, we are interested in assessing the performances of T opt and T orth in transform coding.

Results included in this paper show coding performances on two synthetic data sets and a natural image data set. The synthetic data sets are used to show that T opt and T orth can outperform the KLT when used for high-rate transform coding of non-Gaussian signals. As for the second data set, it consists of well-known grayscale natural images which will be used to evaluate the performances of T opt and T orth in medium-to-low bit rate image compression.

Examples of synthetic signals efficiently compressed with

GCGsup and OrthICA.

The tested signals

The first synthetic data set consists of 2 16 samples of a bidimensional (N = 2) random vector X obtained as follows. First, we produce 2 16 samples of a white vector S = [S 1 , S 2 ] T whose the i-th component is the standardized random variable associated to

S i = Sign(Z i ).|Z i | α , where [Z 1 , Z 2 ]
T is a standardized white Gaussian random vector. The exponent α is an arbitrary positive real number.

When α > 1 (resp. α < 1), S i is super-(resp. sub-) Gaussian. Then, the vector X is obtained via a mirror symmetry, according to the operation X = BS, with B =

Evaluation methodology

Our objective metric to measure the performance of a transform is the generalized coding gain [START_REF] Huang | Block quantization of correlated Gaussian random variables[END_REF].

Given both N sources and a transform, the estimation of ( 5) requires good estimates of the pdfs of each source as well as each transformed component, which may be very difficult to obtain. In this section, we elaborate on a more "practical" way of evaluating [START_REF] Huang | Block quantization of correlated Gaussian random variables[END_REF] which consists in actually coding the transformed components and measuring both the bit-rate-given here by the first order entropy of the quantized data-and the actual end-to-end distortion.

For each tested signal, the vector X is first linearly transformed by a transform T to produce the vector Y = [Y 1 , . . . , Y N ] T whose components are coded separately. The i-th component Y i is first high-rate quantized with a uniform scalar quantizer of quantization step q i . This gives Y i . The bit-rate R i is then estimated by computing the empirical first order entropy of Y i and the inverse transform is

applied to Y = [ Y 1 , . . . , Y N ] T in order to reconstruct an approximation X = [ X 1 , . . . , X N ] T of X. The distortion is the end-to-end one, D = 1 N E[ X -X 2
] and the total average rate is the empirical mean of the N rates R i (1 ≤ i ≤ N ). The optimal allocation of rates between the transform coefficients results in equal weighted distortions w i D i (see section 2). Moreover, using uniform scalar quantizers, bit allocation amounts to choosing a quantization step q i for each of the N components, and for small q i the distortion D i may be well approximated by q 2 i /12, [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. Therefore, the bit allocation follows a simple rule : let c be a constant, then make all the quantization steps q 1 , . . . , q N such that w i D i = c.

This gives q i = 12c/w i , for i = 1, . . . , N . When the constant c varies (under the assumption of high resolution quantization for each component) we obtain the classical asymptotic curve (distortion versus bit-rate, or equivalently bit-rate versus distortion). In our tests, we consider that the hypothesis of high resolution quantization is valid when for each component Y i , the relative deviation between

the actual distortion E[(Y i -Y i ) 2 ]
(where the expectation is estimated by empirical mean) and q 2 i /12

is not greater than 1%. For a given high bit-rate, the ratio between the end-to-end distortion read on the asymptotic curve obtained using the identity transform and that read on the asymptotic curve obtained using T yields the generalized coding gain of T.

Results

The first set of synthetic data was designed so that the KLT does nothing on it. Indeed, the random vector X being white, the generalized coding gain G of the KLT is equal to 0 dB. However, the components of X are not independent, as can be seen on the second row in Fig. 2, and any algorithm among GCGsup, OrthICA and icainf gives the same result: the components Y i are independent (see third row in Fig. 2), and the generalized coding gain G is the same. The generalized coding gain associated with using the transform returned by GCGsup is presented in Tab. 1 for four different values of α. Except for α = 1, i.e., when the components of X are Gaussian, the generalized coding gain of T opt is always greater than zero and becomes more and more important as |1 -α| increases, i.e., as the components of X depart from Gaussianity.

The second set of synthetic data has recently drawn attention of researchers in transform coding.

In [START_REF] Effros | Suboptimality of the Karhunen-Loève transform for transform coding[END_REF], the authors have investigated the strengths and limitations of the KLT when it is used for transform coding the vector X = BU. While the KLT for X is not unique, practical implementations of the KLT give a matrix close to the identity. As a consequence, the diamond distribution of X remains practically unchanged after applying the KLT on X (see Fig. 3(c)). Yet X and U are not equally good sources for scalar quantization as was already pointed out earlier [START_REF] Feng | On the rate-distortion performance and computational efficiency of the Karhunen-Loève transform for lossy data compression[END_REF]. This is confirmed in our experiments since any algorithm among GCGsup, OrthICA and icainf transforms the diamond distribution back into the square distribution of U as can be seen in Fig. 3(d). Furthermore, Tab.

2 shows that the generalized coding gain G of the transform returned by GCGsup outperforms that obtained with the KLT.

Application of GCGsup and OrthICA to the compression of natural images.

The next experimental tests aim at assessing the performances of GCGsup and OrthICA when applied to the compression of well-known grayscale natural images (Lenna, Goldhill, Mandrill, Peppers, Zelda and Boat) each of size 512 × 512 pixels and coded using 8 bits per pixel (bpp) (see Fig. 4). The image coder used in our tests has been designed for experimentation and is not intended to outperform current state-of-the-art image coders such as JPEG2000 [START_REF] Marcellin | An overview of JPEG2000[END_REF]. In natural image compression, we are usually interested in performances at medium-to-low bit rates. In our experiments, we chose to investigate the performances of GCGsup and OrthICA at bit rates less or equal than 2 bpp. Note that for bit-rates less than about 1.6 bpp, the transforms returned by GCGsup and OrthICA can no longer be considered as optimal since the high resolution hypothesis is no longer valid3 .

Bases estimation

The modified ICA bases (i.e., the column vectors of T) were estimated according to two learning schemes. The first scheme yields 12 different bases (2 per image): for each test image, the algorithms (the component which is equivalent to the DC component of the KLT was omitted). The kurtosis was normalized so that it is equal to zero in the case of Gaussian samples. The average of this average kurtosis was computed over the set of test images. The result for each transform is displayed over the label "combined" in Fig. 6. Results show that the kurtosis obtained with the modified ICA bases are greater than those obtained with the KLT and 2-D DCT, regardless of the transform and image considered. Note that the ICA bases exhibit even higher kurtosis, yet yielding the worst generalized coding gains among all tested transforms. The explanation for this is straightforward: ICA bases were estimated using the ICA algorithm icainf which aims at minimizing the mutual information between the components, i.e., only the first term in eq. ( 8). This is done without putting any orthogonality constraint on the basis vectors resulting in a transform which is far from being orthogonal (see Tab. 3).

This result highlights the importance of trading-off between independence and orthogonality-when the distortion is measured using mean squared error. This is well illustrated by eq. ( 8) (a discussion about orthogonality and independence can also be found in [START_REF] Goyal | Theoretical foudations of transform coding[END_REF]).

Medium-to-low bit rates

The image coder used in our experiment is a transform coder originally developed by Davis 4 . It is very modular and allows for simple replacements of individual components (quantizer, entropy coder, transform). It was modified so that it ressembles a JPEG-like coder. The image to be coded is first "tiled" into blocks of eight by eight pixels each, then each tile is represented into a new basis using one of the tested transforms. The bases obtained using the first learning scheme are transmitted with the image since they are data-dependent bases. As for the bases estimated using the second learning scheme, they are not transmitted with the image. Quantization steps are chosen to minimize the end-to-end distortion (2) subject to bit rate constraint. The bit allocation procedure is based on integer programming algorithms described in [START_REF] Shoham | Efficient bit allocation for an arbitrary set of quantizers[END_REF] which provide optimal or near-optimal allocations for the quantizers included here. Entropy coding of the quantizer output is carried out by an adaptive arithmetic coder.

Compression results

We now compare the compression performances of the previously mentioned transforms: KLT, T opt , T orth , KLT , T opt , T orth and 2-D DCT. The results are displayed in Tab. 5, in which we present the peak signal to noise ratio (PSNR) as a function of bit-rate for each test image. Whatever the image, the plots have these common characteristics:

1) The transform codes based on KLT , T opt , T orth and the 2-D DCT perform better than those based on the KLT, T opt and T orth , regardless of the bit-rate. The poor coding performances of the KLT, T opt and T orth are mainly due to the coding penalty resulting from coding the basis vectors (11 bits were allocated on average to each matrix coefficient resulting in a coding precision of 10 -3 ).

In our tests, we observed that the coding penalty is about 1 dB at 2 bpp and 1 bpp, about 2 dB at 0.5 bpp and can reach up to 5 dB at 0.25 bpp. The rate-distortion curves obtained with the KLT, T opt and T orth without coding the basis vectors (not plotted here for the sake of clarity) were close to those obtained with respectively KLT , T opt and T orth (no meaningful performance difference could be observed).

2) At bit-rates greater than 1 bpp, all the rate-distortion curves are lines of slope about 5 ∼ 6 dB per bit. We observed that the high-resolution hypothesis is verified for bit-rates greater than about 1.6 bpp (in this case, the slope of the curve is ∼ 6 dB per bits).

3) Whatever the bit-rate, no meaningful performance difference can be observed between the classadapted transform codes based on T opt and those based on T orth . For medium-to-high bit rates (i.e., greater than about 1.6 bpp), this result can be more or less predicted by looking at the values of the corresponding generalized coding gains (see Tab. 4). For bit-rates greater than 1 bpp, the performance difference between either T opt or T orth and KLT is equal to about the difference between their corresponding generalized coding gains suggesting that the performance of our entropy coder is close to that of a perfect first-order entropy coder. In the low bit-rate region (< 1 bpp), the performance difference tends to become smaller. As for the 2-D DCT, its performance is comparable with that of any of the class-adapted modified ICA bases. Thus, our approach has made it possible to learn two bases which are competitive with the well-known 2-D DCT basis according to the standard PSNR measure.

Tab. 5 compares the performance in compression of the 2-D DCT, the class-adapted KLT and both class-adapted modified ICA bases. The corresponding values of PSNR and achieved bit-rates (the target bit-rate was set equal to 0.5 bpp) are shown. Although the KLT -coded images have worse PSNR than the others, no difference in term of visual quality can be seen. However, when looking further into the details of the reconstructed images Boat in Fig. 7, some meaningful visual quality difference can be seen. Fig. 7) show the images coded with T orth . Black arrows point towards details which are not present or blurred on the corresponding images coded with the 2-D DCT. These details represent lines (e.g., some ropes in the case of the image Boat) which are well preserved with T orth . These results suggest that the class-adapted modified ICA bases are better suited to coding fine details such as lines and edges compared to the 2-D DCT. This is not quite surprising given the more pronounced edge-like nature of the modified ICA bases (see Fig. 5).

Conclusion

This paper addresses the problem of finding optimal 1-D linear transforms in transform coding without

Acknowledgment

The authors are grateful to Pierre Duhamel for his very helpful comments as well as Jacques Weidig and Jean-Louis Gutzwiller who took part in the simulation work.

A Appendix

A.1 Algorithm GCGsup

Given that I(Y 1 , . . . , Y N ) = i h(Y i ) -h(Y) and h(Y) = h(X) + log 2 | det T|, and since the term h(X)

does not depend on T, minimizing the contrast (8) is the same as minimizing

C(T) = C O (T) + C ICA (T)
where

C ICA (T) = N i=1 h(Y i ) -log 2 | det T|.
Using the results of [START_REF] Pham | Entropy of a variable slightly contaminated with another[END_REF] it can be seen that the Taylor expansion of C ICA (T + ET) up to second order may be approximated as follows

C ICA (T + ET) = C ICA (T) + 1≤i =j≤N E[ψ Y i (Y i )Y j ]E ij + 1 2 1≤i =j≤N {E[ψ 2 Y i (Y i )] E[Y 2 j ]E 2 ij + E ij E ji } + • • • , (10)
where the function ψ Y i is equal to the derivative of -log 2 p(y i )p(y i ) denoting the Y i pdf -and is known as the score function. This approximation concerns only the second order terms in the expansion, but not the first order terms. It relies essentially on the assumption of independent transform 20 coefficients, which may not be valid if the solution of the ICA problem is far from the solution that minimizes the contrast [START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF]. But it is quite useful since it leads to a decoupling in the quadratic form of the expansion. Let M = T -T T -1 . One may verify that the Taylor expansion of C O (T + ET) with respect to E and around E = 0, up to second order, is given by

C O (T + ET) = C O (T) -1≤i =j≤N M ji M ii E ji - 1 2 1≤i =j≤N E ij E ji -2 M ji M jj E ii E ji + N i=1 N j=1, j =i N k=1, k =i M jk 2M ii - M ij M ik M 2 ii E ji E ki + M kj M kk E ji E ik + • • • .
The quadratic form associated with the above expansion is quite involved and is not positive.

One possible approximation consists in neglecting the non diagonal elements of M, which amounts to assuming that the optimal linear transform is close to an orthogonal transform. Under this hypothesis, one may verify that

C O (T + ET) ≈ C O (T) - 1≤i =j≤N M ji M ii E ji + 1 2 1≤i =j≤N M jj M ii E 2 ji + E ji E ij + • • • . (11)
The quadratic form associated with the above expansion is now positive, but not positive definite.

However, this is sufficient for the matrix associated with the quadratic form of the Taylor expansion of C(T) to be positive definite, which ensures the stability of the iterative algorithm. Finally, by adding equation [START_REF] Bell | The 'independent components' of natural scenes are edge filters[END_REF] and approximation [START_REF] Puga | An experiment on comparing PCA and ICA in classical transform image coding[END_REF] we obtain an approximation of C(T + ET) up to second order of E ij and the iteration consists explicitly of solving the linear equations

    E[ψ 2 Y i (Y i )]E[Y 2 j ] + M ii M jj 2 2 E[ψ 2 Y j (Y j )]E[Y 2 i ] + M jj M ii         E ij E ji     =     M ij M jj -E[ψ Y i (Y i )Y j ] M ji M ii -E[ψ Y j (Y j )Y i ]     .
The indeterminate diagonal terms E ii are arbitrarily fixed to zero. Then the estimator T is left multiplied by I + E in order to update it. In this expression, the probability density functions being unknown, the score function ψ Y i (y i ) is replaced by an estimation (see [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF]) and the expectations are estimated by empirical means.

A.2 Algorithm OrthICA

In this section, we propose to find the orthogonal transform that minimizes the contrast [START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF]. Since the second term of (8) vanishes for any othogonal matrix T, this amounts to finding the orthogonal transform which minimizes the first term of (8), or equivalently, which minimizes C ICA (T). If the matrix T is orthogonal, so is T + ET, providing that I + E be orthogonal. This last condition will be satisfied up to second order if E is anti-symmetric, since (I + E) T (I + E) = I + E T E differs from the identity only by second order terms. Let E be anti-symmetric. The Taylor expansion of

C ICA (T + ET) becomes C ICA (T + ET) = C ICA (T) + 1≤i<j≤N E[ψ Y i (Y i )Y j ] -E[ψ Y j (Y j )Y i ] E ij + 1 2 1≤i<j≤N E 2 ij E[ψ 2 S i (Y i )] E[Y 2 j ] + E[ψ 2 Y j (Y j )] E[Y 2 i ] -2 + • • •
, and the minimization of the second term in the above expansion yields

E ij = E[ψ Y j (Y j )Y i ] -E[ψ Y i (Y i )Y j ] E[ψ 2 Y i (Y i )] E[Y 2 j ] + E[ψ 2 Y j (Y j )] E[Y 2 i ] -2 . (12)
Actually, T + ET is not a true orthogonal transform. This may be overcome by replacing T + ET with e E T = (I + E + E 2 /2! + • • • )T, which is an orthogonal matrix differing from T + ET only by second order terms.

Entropy coding

) ( ˆm x N ) ( ˆ2 m x ) ( ˆ1 m x ) ( ˆm y N ) ( ˆ2 m y ) ( ˆ1 m y ) ( 1 m y ) ( 2 m y ) (m y N U ) ( x m ˆ ) ( 1 m x ) ( 2 m x ) (m x N T Q 1 Q 2 Q N ) (m x E 1 E 2 E N Analysis transform Synthesis transform Quantization ) ( x n ) ( x ˆn

3. 1 Criteria

 1 The criterion[START_REF] Jana | Optimality of KLT for high-rate transform coding of Gaussian vectorscale mixtures: application to reconstruction, estimation and classification[END_REF] may be decomposed into C(T) = C ICA (T) + C O (T), where C ICA (T) = I(Y 1 ; . . . ; Y N ) corresponds to the mutual information criterion in ICA, and

  matrix. Samples of S and X are depicted respectively on the first row and the second row in Fig. (2) for four different values of α. The second synthetic data set consists of 2 16 samples of a bidimensional random vector X obtained as follows: X = BU, where B is given above and U = [U 1 , U 2 ] T with U 1 and U 2 independent and identically distributed uniform random variables on [-1, 1]. Samples of U and X are depicted respectively in Fig. 3(a) and Fig. 3(b).

GCGsup

  and OrthICA were applied to a training set consisting of 4 096 non overlapping image blocks each of size 8 × 8 pixels extracted from the test image. The first and second columns of Fig.5displays the estimated modified ICA bases as well as the practically achieved KLT bases obtained from the test images Boat and Peppers, respectively. Also displayed for comparison are the bases obtained using the ICA algorithm icainf. As for the second scheme, it yields only 2 different bases. The modified ICA bases were learned from one training set consisting of 12 288 non overlapping image blocks each of size 8 × 8 pixels extracted from three test images (Lenna, Goldhill and Boat). The third column of Fig.5displays the estimated modified ICA bases as well as the KLT basis (denoted KLT ). The transform whose column vectors were estimated using the algorithm OrthICA (resp. GCGsup) is denoted T orth (resp. T opt ). The ICA basis obtained with the icainf algorithm for this training set is also displayed for comparison and denoted T ICA .Examining Fig.5closely reveals that the features found with GCGsup and OrthICA are much more localized in space than the checkerboardlike basis vectors obtained with the KLT. Notice also the more transform coefficients of the KLT and 2-D DCT, which tend to a normal distribution. Fig.6displays for each transform and each image the average kurtosis computed over 63 transformed components

  the classical assumption of Gaussianity. This paper emphasizes a new point of view in variable rate transform coding by showing, under the high resolution hypothesis, that the problem of finding the optimal 1-D linear transform may be recast as a modified independent component analysis (ICA) problem. Two new modified ICA algorithms, called GCGsup and OrthICA, are introduced for computing the optimal 1-D linear transform and the optimal 1-D orthogonal transform, respectively.Experimental results included in this paper show coding performances of GCGsup and OrthICA on two synthetic data sets and a natural image data set. Experiments carried out with the synthetic data sets show that the new transforms can outperform the KLT when used for high-rate transform coding of non-Gaussian signals. When applied to the compression of some well-known natural images, GCGsup and OrthICA have proved 1) to be comparable to the classical 2-D DCT according to the PSNR measure and 2) to yield better visual image quality (better preservation of lines and edges) than the 2-D DCT.

Figure 1 :

 1 Figure 1: Compression system including a linear transform, a quantization and an entropy coding stage.

  (a) S; α = 1 (b) S; α = 1.5 (c) S; α = 2 (d) S; α = 2.5(e) X (f) X (g) X (h) X (i) Y, Topt (j) Y, Topt (k) Y, Topt (l) Y, Topt

Figure 2 :

 2 Figure 2: Samples of S, X and Y (obtained after applying T opt on each sample of X) for four different values of α.

Figure 3 :Table 1 :Table 2 :

 312 Figure 3: Samples of U, X and Y (obtained after applying the KLT and T opt on X).

Figure 4 :

 4 Figure 4: Test images.

Figure 5 :

 5 Figure 5: KLT, T orth , T opt and T ICA basis vectors obtained from boat and peppers (on first and second column, respectively) and KLT , T orth , T opt and T ICA basis vectors obtained from the image class training set. 27

Figure 6 :

 6 Figure 6: Average kurtosis computed over 63 transformed components (the component which is equivalent to the DC component of the KLT was omitted).

  29

  T orth (f) T orth (g) T orth (h) T orth

Figure 7 :

 7 Figure 7: Zoom on parts of image Boat coded at about 0.5 bpp. Black arrows point towards details which are not present or blurred on the corresponding image coded with the 2-D DCT.

Table 3 :

 3 Distance to orthogonality (in bpp) of T opt , T opt , T ICA and T ICA for each test image.

Table 4 :

 4 Generalized coding gain (in dB) of the KLT, T orth , T opt , KLT , T orth , T opt , T ICA , T ICA and 2-D DCT for each test image. Last column yields the average generalized coding gain of each transform computed over the set of test images.

		Lenna Goldhill Mandrill Peppers Zelda Boat Average
	KLT	18.13	15.22	6.99	17.23	19.91 15.35	15.47
	T orth	18.58	15.62	7.42	17.89	20.12 15.96	15.93
	T opt	18.60	15.68	7.49	17.84	20.12 15.94	15.94
	KLT	17.76	15.09	6.98	17.08	19.26 14.79	15.16
	T orth	18.34	15.40	7.27	17.63	19.89 15.74	15.73
	T opt	18.31	15.46	7.28	17.70	19.84 15.72	15.71
	T ICA	6.88	11.63	5.52	11.47	5.36 11.62	8.74
	T ICA	16.50	13.75	5.51	16.03	17.83 13.95	13.92
	2-D DCT 18.25	15.37	7.17	17.50	19.87 15.61	15.62

Table 5 :

 5 PSNR (dB) versus target bit-rate for each test image and each transform.

As in Section

4.1.2, we consider that the hypothesis of high resolution quantization is valid when for each component Yi, the relative deviation between the actual distortion E[(Yib Yi) 2 ] (where the expectation is estimated by empirical mean) and q 2 i /12 is not greater than 1%.

http://www.geoffdavis.net/dartmouth/wavelet/wavelet.html

This work was supported in part by the French Ministry of Higher Education and Research, with the "Action Concertée Incitative Masse de Données" ACI 2 M project, and * by the Lorraine Region.

Other similar experiments with ICA [START_REF] Hyvärinen | Independent Component Analysis[END_REF] have produced comparable results.

The computation of the pseudo-distance to orthogonality (see eq. ( 9)) reveals that, unlike the bases obtained with the ICA algorithm icainf, those estimated with GCGsup are quasi-orthogonal as can be seen in Tab. 3. As mentioned earlier, it is not surprising that orthogonal transforms are favored over other linear transforms since we measure the distortion as mean squared error and orthogonal transforms are energy-preserving.

Medium-to-high bit rates

Tab. 4 shows estimations of the generalized coding gain for each tested transform and each test image.

The average generalized coding gain computed over the set of test images is also given. The estimation method used is the same as that described in section 4.1.2. Looking at the average values of the generalized coding gain reveals that, whatever the learning scheme, the modified ICA transforms perform best followed respectively by the 2-D DCT, the KLT and, far behind, the ICA transform. The coding gain of any of the modified ICA transforms relative to the 2-D DCT is about 0.3 dB (resp. 0.1 dB) when the first (resp. second) learning scheme is used suggesting that a transform-based image coder could benefit from using any of the modified ICA transforms. However, these results should be taken with care since they are meaningful only under the high resolution hypothesis (i.e., in general, at medium-to-high bit rates). Yet, in image compression, we are most interested in performances at medium-to-low bit rates, where the high resolution hypothesis is no longer valid. Additional precaution should be taken for the bases estimated according to the first learning scheme. In this case, the basis vectors need to be transmitted (coded) with the image resulting in an extra coding cost which has not been incorporated in the method used for estimating the generalized coding gain.

Another motivation for using the modified ICA bases in compression lies in the transform coefficients distribution which tends to be heavy-tailed, more suited for quantization and entropy coding than the