Ahmed Frikha

Patrice Cartraud

Fabien Treyssède
email: fabien.treyssede@ifsttar.fr

Mechanical modeling of helical structures accounting for translational invariance. Part 1: Static behavior

Keywords: Homogenization, Helical coordinates, Finite element method, Helical springs, Seven-wire strands, Axial load

The purpose of this paper is to investigate the static behavior of helical structures under axial loads. Taking into account their translational invariance, the homogenization theory is applied. This approach, based on asymptotic expansion, gives the first-order approximation of the 3D elasticity problem from the solution of a 2D microscopic problem posed on the cross-section and a 1D macroscopic problem, which turns out to be a Navier-Bernoulli-Saint-Venant beam problem. By contrast with earlier references in which a reduced 3D model was built on a slice of the helical structure, the contribution of this paper is to propose a 2D microscopic model. Homogenization is first applied to helical single wire structures, i.e. helical springs. Next, axial elastic properties of a sevenwire strand are computed. The approach is validated through comparison with reference results: analytical solution for helical single wire structures and 3D detailed finite element solution for seven-wire strands.

Introduction

Helical structures are widely used in mechanical and civil engineering applications. These structures are usually subjected to large loads which can lead to the material degradation and cracks associated with corrosion and mechanical fatigue. This threatens the structural strength. In this framework, nondestructive testing is a crucial tool for detection, localisation and measurement of material discontinuities. The choice of the appropriate technique depends on dimensions and accessibility of the structure. Particularly, ultrasonics allow to control large components, such as plates and tubes, by analyzing their elastic guided waves. The purpose of this study, which is composed of two parts, is to develop a numerical model for the analysis of the elastic wave propagation phenomenon in prestressed helical structures. This problem requires the computation of the static prestress state. Therefore, a first model will be developed in Part 1 of this paper, to compute this static state. Taking into account this prestress state, a second model will be developed in Part 2, in order to analyze the wave propagation in these prestressed structures. The goal of this first part of this paper is thus to develop an approach that allows the computation of the prestress state in helical structures subjected to axial load.

Numerous works have been devoted to the modeling of the static behavior of helical structures as springs and multi-wire cables under axial loads. For helical springs, an analytical model was proposed among others in [START_REF] Ancker | Pitch and curvature corrections for helical springs[END_REF] and [START_REF] Wahl | Mechanical Springs. Second Edition[END_REF] considering the spring as an Euler-Bernoulli beam with pitch and curvature corrections. Numerical approaches describing the static behavior of helical springs have been also developed. Among these works, a finite element model of half of a spring slice has been proposed in [START_REF] Jiang | A novel finite element model for helical springs[END_REF].

The static behavior of seven-wire strands has been widely studied in literature. Various analytical models based on different assumptions have been proposed, such as the model of [START_REF] Costello | Theory of Wire Rope[END_REF] which is one of the most popular. These models are reviewed in [START_REF] Jolicoeur | A numerical comparison of current mathematical models of twisted wire cables under axisymmetric loads[END_REF] and compared in [START_REF] Jolicoeur | A numerical comparison of current mathematical models of twisted wire cables under axisymmetric loads[END_REF] and [START_REF] Ghoreishi | Validity and limitations of linear analytical models for steel wire strands under axial loading, using a 3d fe model[END_REF]. Besides, numerical models relying on the finite element method were developed. Some of them are based on beam elements [START_REF] Durville | Modélisation du comportement mécanique de câbles métalliques[END_REF]; [START_REF] Nawrocki | A finite element model for simple straight wire rope strands[END_REF]; [START_REF] Páczelt | Nonlinear contact-theory for analysis of wire rope strand using high-order approximation in the fem[END_REF]), see also [START_REF] Nemov | Generalized stiffness coefficients for iter superconducting cables, direct fe modeling and initial configuration[END_REF] and [START_REF] Bajas | Numerical simulation of the mechanical behavior of iter cable-in-conduit conductors[END_REF] in which ITER superconducting cables composed of a large number of strands are studied. But most of the time, 3D models are used, see e.g. [START_REF] Boso | Homogenisation methods for the thermo-mechanical analysis of nb3sn strand[END_REF], [START_REF] Ghoreishi | Validity and limitations of linear analytical models for steel wire strands under axial loading, using a 3d fe model[END_REF], [START_REF] İmrak | On the problem of wire rope model generation with axial loading[END_REF], [START_REF] Nemov | Generalized stiffness coefficients for iter superconducting cables, direct fe modeling and initial configuration[END_REF], Stanova et al. (2011a,b), [START_REF] Erdönmez | New approaches for model generation and analysis for wire rope[END_REF]. In order to obtain a good representation of the geometry as well as the displacement solution, which may involve bending phenomena, quadratic elements are employed. This leads to models which can be computationally expensive, when the model axial length is about the pitch length. Therefore, as soon as the loading fulfills helical symmetry, one can take benefit of this property to reduce the model size. This has been achieved in [START_REF] Jiang | A concise finite element model for simple straight wire rope strand[END_REF][START_REF] Jiang | Statically indeterminate contacts in axially loaded wire strand[END_REF] in which the computational domain is restricted to a basic sector of a helical slice. Helical symmetry may also be accounted for within the framework of homogenization theory. This has been proposed first in [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF] using axial periodicity, and then improved in [START_REF] Messager | Homogenization of helical beam-like structures: application to single-walled carbon nanotubes[END_REF], in which helical symmetry enables to consider one slice of a strand. The derivation of the slice model is different in [START_REF] Jiang | A concise finite element model for simple straight wire rope strand[END_REF][START_REF] Jiang | Statically indeterminate contacts in axially loaded wire strand[END_REF] and [START_REF] Messager | Homogenization of helical beam-like structures: application to single-walled carbon nanotubes[END_REF]. However, in both cases, helical symmetry yields displacement constraints between the two faces of the slice, with a loading under the form of an axial strain and a twist rate. This work further advances [START_REF] Cartraud | Computational homogenization of periodic beam-like structures[END_REF] and [START_REF] Messager | Homogenization of helical beam-like structures: application to single-walled carbon nanotubes[END_REF], taking advantage of the translational invariance. Helical symmetry can be actually considered more efficiently. Thus the model can be reduced to a 2D one, i.e. a cross-section model. This requires to formulate the homogenization theory in a twisted coordinate system. This technique then allows the computation of the static prestressed state of helical structures (single wire and multi-wire) from the solution of a 2D problem. Let us mention that an advanced analytical 2D model has been recently proposed in [START_REF] Argatov | Response of a wire rope strand to axial and torsional loads: Asymptotic modeling of the effect of interwire contact deformations[END_REF].

This model takes into account Poisson's effect, contact deformation and allows to obtain the overall strand stiffness as well as local contact stresses. In this reference, plane strain was assumed to formulate the 2D problem while in the present work helical symmetry is used.

The method developed in this paper is restricted to multi-wire helical struc-tures composed of a stack of helical wires wrapped with the same twisting rate around a straight axis. As explained in Section 3, this excludes the case of double helical structures (such as independent wire rope core for instance) and cross-lay strands. This paper is organized as follows. First, in Section 2, the curvilinear coordinate system is introduced. Then in Section 3 the translational invariance is defined, which is a necessary condition for the helical homogenization approach.

Based on the asymptotic expansion method and exploiting the translational invariance property, the homogenization procedure is presented in Section 4.

Its finite element solution is detailed in Section 5. The helical homogenization approach is validated for helical single wire and seven-wire structures by comparison with analytical or numerical models in Section 6.

Curvilinear coordinate system

A helical structure is considered (see Fig. 1). Let (e X , e Y , e Z ) its Cartesian orthonormal basis. The helix centreline is defined by its helix radius R in the Cartesian plane (e X , e Y ) and the length of one helix pitch along the Z-axis denoted by L. This helix centerline can be described by the following position vector:

r(s) = R cos( 2π l s + θ)e X + R sin( 2π l s + θ)e Y + L l se Z , (1)
where l = √ L 2 + 4π 2 R 2 is the curvilinear length of one helix pitch and θ is the helix phase angle in the Z = 0 plane. For a seven-wire strand, θ is equal to (N -1)π/3, where N = 1, .., 6 refers to the number of the helical wire. θ is equal to zero for a single wire helical structure. The helix lay angle Φ is defined by tanΦ = 2πR/L. A complete helix is described by the parameter s varying from 0 to l.

Serret-Frenet basis

A Serret-Frenet basis (e n , e b , e t ) associated to the helix can be defined (see e.g. [START_REF] Gray | Modern Differential Geometry of Curves and Surfaces with Mathmatica[END_REF], where the unit vectors e n , e b , e t are given by e t = dr/ds, de n /ds = τ e b -κe t and de b /ds = -τ e n . For helical curves, the curvature κ = 4π 2 R/l 2 and the torsion τ = 2πL/l 2 are constant. In the Cartesian basis, e n , e b and e t are expressed by:

e n = -cos( 2π l s + θ)e X -sin( 2π l s + θ)e Y , e b = L l sin( 2π l s + θ)e X - L l cos( 2π l s + θ)e Y + 2π l Re Z , e t = - 2πR l sin( 2π l s + θ)e X + 2πR l cos( 2π l s + θ)e Y + L l e Z .
(2)

The normal vector e n remains parallel to the (e X , e Y ) plane while e b and e t move in the three directions of the Cartesian basis as s and θ vary.

Twisted basis

A special case of the Serret-Frenet basis denoted by (e x , e y , e Z ) corresponding to κ = 0 and τ = 2π/L can be considered. It corresponds to a twisted coordinate system along the Z-axis (s ≡ Z) with axial periodicity L. The unit vectors e x and e y rotate around the Z-axis and remain parallel to the (e X , e Y ) plane (see Fig. 1). In the Cartesian basis, e x and e y are expressed as:

e x = -cos( 2π L Z + θ)e X -sin( 2π L Z + θ)e Y , e y = sin( 2π L Z + θ)e X -cos( 2π L Z + θ)e Y .
(3)

It should also be noted that this twisted coordinate system coincides with the one proposed in [START_REF] Onipede | Propagating waves and end modes in pretwisted beams[END_REF], [START_REF] Nicolet | Modeling of twisted optical waveguides with edge elements[END_REF][START_REF] Nicolet | Finite element analysis of helicoidal waveguides[END_REF] for the analysis of twisted and helical structures.

Covariant and contravariant bases

Differential operators can not be expressed directly in the Serret-Frenet or twisted bases. They have first to be expressed in the covariant and contravariant bases. The reader can find an in-depth treatment of curvilinear coordinate systems in [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF], [START_REF] Synge | Tensor Calculus[END_REF], [START_REF] Wempner | Mechanics of Solids with Applications to Thin Bodies[END_REF] for instance. From the twisted basis (e x , e y , e Z ), a new coordinate system (x, y, Z) is built, for which any position vector can be expressed as:

X(x, y, Z) = xe x (Z) + ye y (Z) + Ze Z . (4)
The covariant basis (g 1 , g 2 , g 3 ) is obtained from the position vector by (g 1 , g 2 , g 3 ) =

(∂X/∂x, ∂X/∂y, ∂X/∂Z), which yields in the twisted basis:

g 1 = e x (Z) , g 2 = e y (Z) , g 3 = -τ ye x (Z) + τ xe y (Z) + e Z .
(5)

Note that the covariant basis is not orthogonal.

The covariant metric tensor, defined by g mn = g m • g n , is then given by:

g =      1 0 -τ y 0 1 τ x -τ y τ x τ 2 (x 2 + y 2 ) + 1      . ( 6
)
The covariant basis gives rise to the contravariant one (g 1 , g 2 , g 3 ), defined from g i • g j = δ j i . Superscripts and subscripts refer to the covariant and contravariant vectors, respectively. g 1 , g 2 and g 3 are expressed in the twisted basis as:

g 1 = e x (Z) + τ ye Z , g 2 = e y (Z) -τ xe Z , g 3 = e Z . (7)
The Christoffel symbol of the second kind Γ k ij , defined by Γ k ij = g i,j • g k , can be calculated from the covariant and contravariant bases, which leads to:

Γ k 11 = Γ k 12 = Γ k 21 = Γ k 22 = 0, Γ 1 13 = Γ 1 31 = 0, Γ 1 23 = Γ 1 32 = -τ, Γ 1 33 = -τ 2 x, Γ 2 23 = Γ 2 32 = 0, Γ 2 33 = -τ 2 y, Γ 2 13 = Γ 2 31 = τ, Γ 3 13 = Γ 3 31 = Γ 3 23 = Γ 3 32 = Γ 3 33 = 0. ( 8
)
It is noteworthy that the coefficients Γ k ij do not depend on the axial variable Z. As shown in the next section, this is a necessary condition for translational invariance.

Strain tensor

The strain tensor is now rewritten in the curvilinear coordinate system. In the contravariant basis, the strain-displacement relation is [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF]):

ǫ = ǫ ij g i ⊗ g j , ǫ ij = 1 2 (u i,j + u j,i ) -Γ k ij u k , (9)
where the u i 's denote the displacement covariant components.

Using the relation ( 7) between the contravariant and the twisted bases, the strain vector can then be expressed in the twisted basis as follows:

{ǫ} = (L xy + L Z ∂ ∂Z ){u}, L xy =               ∂/∂x 0 0 0 ∂/∂y 0 0 0 Λ ∂/∂y ∂/∂x 0 Λ -τ ∂/∂x τ Λ ∂/∂y               , L Z =               0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0               , (10)
where Λ = τ (y∂/∂x -x∂/∂y). The column vectors {u} = [u x u y u Z ] T and {ǫ} = [ǫ xx ǫ yy ǫ ZZ 2ǫ xy 2ǫ xZ 2ǫ yZ ] T are the displacement vector and the strain vector respectively, both written in the orthonormal twisted basis (e x , e y , e Z ).

Constitutive law

For an isotropic material, the elasticity tensor is given in the covariant basis by [START_REF] Chapelle | The Finite Element Analysis of Shells-Fundamentals[END_REF]):

C = C ijkl g i ⊗ g j ⊗ g k ⊗ g l , C ijkl = νE (1 + ν)(1 -2ν) g ij g kl + E 2(1 + ν) (g ik g jl + g il g jk ), (11)
where E and ν are the Young modulus and the Poisson's ratio, respectively.

Using the relation between the covariant and the twisted bases and after simplifications, it can be checked that the elasticity tensor components in the twisted basis are given by:

C αβδγ = νE (1 + ν)(1 -2ν) δ αβ δ δγ + E 2(1 + ν) (δ αδ δ βγ + δ αγ δ βδ ), (12)
where greek subscripts {α, β, γ, δ} denote components x, y, Z in the twisted basis. The above expression coincides with the one obtained in the Cartesian basis, as the twisted basis is orthonormal.

Translational invariance

Translational invariance is a key property for applying the homogenization theory. For cylindrical structures, translational invariance means that both the cross-section and the material properties do not vary along the axis. For curved structures, there is another condition which states that the differential operator coefficients must not depend on the axial variable [START_REF] Treyssède | Mode propagation in curved waveguides and scattering by inhomogeneities: application to the elastodynamics of helical structures[END_REF]).

As a consequence, for helical structures, the translational invariance requires the following three conditions [START_REF] Treyssède | Elastic waves in helical waveguides[END_REF], Treyssède and Laguerre ( 2010)):

1. The material properties do not vary along the Z-axis in the twisted coordinate system;

2. The coefficients of the differential operators (gradient, divergence, Laplacian, ...) are independent on the axial variable Z;

3. The cross-section does not vary along the Z-axis in the twisted coordinate system.

Throughout this work, the material is assumed to be homogeneous and isotropic. In this case, the first condition is verified. To satisfy the second condition, it is sufficient to prove that the Christoffel symbols do not depend on the axial variable Z, which has been verified in the last section (see Eq. 8).

Thus it remains only to verify the third condition.

Let us consider a helical single wire structure. The cross-section shape in the (e X , e Y ) plane at the axial position Z 1 is similar to that given at the position Z 2 : there only exists a rotation of angle 2π(Z 2 -Z 1 )/L around the Z-axis between these two cross-section shapes. Moreover, because the twisted basis plane (e x , e y ) also rotates around Z, the cross-section indeed remains fixed in this plane. Therefore, the translational invariance is checked for helical single wire structures. Fig. 2 shows the cross-section of four helical single wires with R = 2a and different helical angles in the (e X , e Y ) plane. a is the radius of the circular cross-section (the cross-section being circular in the plane normal to the helical curve). Note that for small angle Φ, the cross-section shape in this plane is nearly circular because the structure is close to a cylinder (Fig. 2a). However the cross-section shape deviates from the circular one as Φ increases.

Let us now consider multi-wire helical structures. They are composed of a stack of helical wires, wrapped around a straight wire. A seven-wire strand is a special case of helical multi-wire structures containing one layer of six helical wires wrapped around the central wire. In the twisted basis, a cylindrical structure of axis Z with isotropic material is translationally invariant for any value of the torsion τ (see [START_REF] Treyssède | Investigation of elastic modes propagating in multi-wire helical waveguides[END_REF]). It therefore remains fixed in the Cartesian as in the twisted coordinate systems. The central wire is hence translationally invariant. As shown for helical single wire, the peripheral helical wires, which have the same helix parameters are also translationally invariant in the twisted coordinate system. The geometric invariance is then verified for the seven-wire strand in the twisted coordinate system and the problem is translationally invariant.

Let us briefly examine more complex structures. In multi-layer wire ropes, more than one layer of helical wires is present. Translational invariance in such structures is still satisfied if the torsion of each wire remains identical. This implies that translational invariance is not fulfilled in case of cross-lay strands because the torsion can be positive or negative. This loss of invariance is obvious if one thinks of contact discontinuities between two layers of opposite torsion.

Contact discontinuities also necessarily occur in double helical structures, composed of one central strand wrapped by several peripheral strands. Such double helical structures, sometimes referred to as IWRC (independent wire rope core), hence cannot fulfill translational invariance.

To conclude this section, let us define the cross-section boundary in the plane Z = 0. The surface boundary of a helical single wire with circular cross-section is described in the Serret-Frenet basis by the following position vector:

X(x, y, s) = r(s) + a cos te n (s) + a sin te b (s), ( 13
)
where t ∈ [0; 2π]. Substituting Eq. ( 2) into Eq. ( 13), the cross-section shape parameterization in the (e X , e Y ) plane is:

               X(t) = (R -a cos t) cos(ηa sin t + θ)+ L l a sin t sin(ηa sin t + θ) Y (t) = (R -a cos t) sin(ηa sin t + θ)+ L l a sin t cos(ηa sin t + θ) , (14)
where η = -4π 2 R/lL. This curve has been used to plot the cross-sections on Fig. 2. It has also been used for the FE mesh generation in Section 6.

Helical homogenization procedure

In this work, helical structures are supposed to be subjected to external loads at its end sections. Moreover, only axial loads (traction and torsion) are considered. Targeted helical structures are helical springs and seven-wire strands.

As explained in introduction, the purpose of this paper is to propose an approach for obtaining the static stress state, which will be used in the second part of this paper as a prestress state, for a wave propagation analysis. This can be achieved efficiently using an homogenization method. This approach, based on the asymptotic expansion method, exploits the translational invariance property. Homogenization splits the initial 3D elasticity problem into 2D problems posed on the cross-section, and a 1D straight beam problem. The overall beam behavior is computed thanks to the solution of the 2D problems. This solution, combined with the solution of the beam problem, provides also the local stress state.

For the present work, let us consider a slender helical structure of axial length H (see Fig. 3), with a cross-section denoted S ε . This structure occupies

the configuration Ω ε = S ε × [0, H]. The boundary of Ω ε is defined by ∂Ω ε = Γ ε ∪Γ ε 0 ∪Γ ε H , with Γ ε 0 = S ε ×{0} and Γ ε H = S ε
×{H} the two end cross-sections of the helical structure and Γ ε the cross-section boundary. This structure exhibits a small parameter ε, corresponding to the inverse of the slenderness ratio, i.e. the ratio between the diameter of the cross-section S ε and the length H. a)

 H    e Z S  H  0  b) H    H   0  Z 1 Z 2 0 Z

The initial problem

The linear elasticity problem consists in finding the fields σ ε , ǫ ε and u ε , solution of:

               ∇ • σ ε = 0 σ ε = C : ǫ ε (u ε ) ǫ ε (u ε ) = ∇ s (u ε ) σ ε • n = 0 on Γ ε , ( 15
)
where C is the elasticity tensor, which is supposed to be constant under the assumption of small displacements. ∇ s (•) and ∇ • (•) denote respectively the symmetric gradient (strain) and divergence operators. The solution must also verify the boundary conditions at the end sections. They are supposed to be under the form of stress data: σ ε • (-e Z ) = t 0 on Γ ε 0 and σ ε • e Z = t H on Γ ε H , where t 0 and t H are the tractions at the end sections located at Z = 0 and Z = H. Moreover t 0 and t H are such that the overall structure equilibrium is fulfilled, which is a necessary condition for problem (15) to have a solution.

For seven-wire strand, the solution must verify Eq. ( 15) on each wire as well as contact equations, on the contact line between the central wire and each helical wire. This raises the problem of contact assumptions. In [START_REF] Ghoreishi | Validity and limitations of linear analytical models for steel wire strands under axial loading, using a 3d fe model[END_REF], stick and slip conditions have been studied for computing the overall behavior. In [START_REF] Gnanavel | Effect of interfacial contact forces in radial contact wire strand[END_REF], an analytical model with frictional contact was developed. Overall stiffness as well as maximum normal contact stresses were calculated from the authors's model and the Costello's model which assumes stick contact. In [START_REF] Argatov | Response of a wire rope strand to axial and torsional loads: Asymptotic modeling of the effect of interwire contact deformations[END_REF] hypothesis of slip contact is made and maximum contact pressures (core-wire and wire-wire) were compared to FE computations performed with frictional contact in [START_REF] Jiang | Statically indeterminate contacts in axially loaded wire strand[END_REF]. All of these previous works have shown that the overall stiffness and contact stresses are very little sensitive to contact conditions. Therefore in this work, for simplicity, the contact is assumed to be stick. This amounts to perfect bonding conditions between wires: u c = u p and (σ

• n) + c + (σ • n) - p = 0
, where the subscripts c and p are related to the central and peripheral wires.

The solution of this problem (15) with boundary conditions and contact equations for multi-wire strand provides the prestress state. As mentioned previously, this problem may be computationally expensive to solve under this form, and the homogenization method aims to simplify it.

Asymptotic expansion method

To our knowledge very few works have been devoted to the asymptotic analysis of helical structures starting from a 3D formulation. We just mention Nicolet et al. (2007) in the framework of electrostatics. Therefore, the approach presented in this paper is based on [START_REF] Buannic | Higher-order asymptotic model for a heterogeneous beam, including corrections due to end effects[END_REF] and Buannic and Cartraud (2001a) developed for axially invariant and periodic beam-like structures respectively. More about asymptotic expansion method for slender structures may be found in some books (Sanchez-Hubert and Sanchez-Palencia (1992); [START_REF] Kalamkarov | Analysis, Design and Optimization of Composite Structures[END_REF]; [START_REF] Trabucho | Mathematical modelling of rods. Handbook of Numerical Analysis IV[END_REF]).

The first step of the method consists in defining a problem equivalent to problem (15), but posed on a fixed domain that does not depend on the small parameter ε. A change of variables is thus introduced which takes into account the structure slenderness, in the twisted coordinates system: (x, y, ζ) = (x, y, εZ). Next, the displacement solution is searched under an asymptotic expansion form:

u(x) = u 0 x (ζ)e x + u 0 y (ζ)e y + εu 1 (x, y, ζ) + ε 2 u 2 (x, y, ζ) + ... (17)
In this expression, the translational invariance is taken into account since the k th -order displacement u k (x, y, ζ) does not depend on the microscopic axial coordinate Z. Moreover, it is usually considered that the 0 th -order displacement has no axial component, which results from the property that for slender structures, the bending stiffness is much lower than axial stiffness. So 0 th -order displacement corresponds to a transverse deflection. Note that a proof of this result may be found in [START_REF] Trabucho | Mathematical modelling of rods. Handbook of Numerical Analysis IV[END_REF] for homogeneous beams, and in [START_REF] Kolpakov | Calculation of the characteristics of thin elastic rods with a periodic structure[END_REF] for beams with periodic structure. As axial loads are considered in this work, and under the assumption that bending is not coupled with tension and torsion, this 0 th -order term vanishes.

Reporting expansion (17) in problem (15) with the use of ( 16), and considering ζ and {x, y} as independent coordinates, one is led to a sequence of problems. On one hand 2D microscopic problems posed on the cross-section S, which will be denoted P m 2D , where m denotes the order of ε in the equilibrium equation. On the other hand a sequence of 1D macroscopic problems will be also obtained, but only the lowest order macroscopic problem will be considered in the following.

Microscopic problems

The lowest order 2D microscopic problem posed on the cross-section S is P 1 2D with the following equations:

               ∇ xy • σ 1 = 0 σ 1 = C : ǫ 1 ǫ 1 = ∇ s xy (u 1 ) σ 1 • n = 0 on ∂S . ( 18
)
It is important to notice that though this problem is 2D, the displacement u 1 has three components. This results from the property than in a matrix form, from Eq. ( 10), one has:

{∇ s xy (u 1 )} = L xy {u 1 } = L xy          u 1 x u 1 y u 1 ζ          . ( 19
)
Problem P 1 2D is well-posed and has a unique solution up to a rigid body motion (Sanchez-Hubert and Sanchez-Palencia (1992); [START_REF] Buannic | Higher-order asymptotic model for a heterogeneous beam, including corrections due to end effects[END_REF]).

The stress solution is obviously equal to zero. The displacement is thus a rigid body motion solution of ∇ s xy (u 1 ) = 0, its expression in the twisted basis is:

u 1 = u 1 ζ (ζ)e Z + ϕ 1 (ζ)[xe y -ye x ], (20)
corresponding to an overall translation u 1 ζ and rotation ϕ 1 around the Z-axis. The solution of problem ( 18) is then given by u 1 with at this step arbitrary

u 1 ζ (ζ) and ϕ 1 (ζ) and ǫ 1 = σ 1 = 0.
The next order microscopic problem P 2 2D involves σ 2 , ǫ 2 and u 2 solution of:

               ∇ xy • σ 2 = 0 σ 2 = C : ǫ 2 ǫ 2 = ∇ s xy (u 2 ) + ∇ s ζ (u 1 ) σ 2 • n = 0 on ∂S . ( 21
)
Note that the displacement vector u 1 , obtained from the solution of the problem P 1 2D , appears through ∇ s ζ (u 1 ) in P 2 2D . From Eq. ( 20), it can be seen that the components of this strain tensor are, under a matrix form:

{∇ s ζ (u 1 )} = 0 0 E E 0 -yE T xE T T , (22)
where E E = ∂u 1 ζ /∂ζ and E T = ∂ϕ 1 /∂ζ and thus can be identified as macroscopic strains, i.e. extension and torsion respectively. Therefore, the other part of the strain ∇ s xy (u 2 ) is a microscopic strain. Thanks to the problem linearity, its solution is a linear function of the macroscopic strains, up to a rigid body displacement which is of the form (20). So one has:

u 2 = χ E (x, y)E E (ζ) + χ T (x, y)E T (ζ)+ u 2 ζ (ζ)e Z + ϕ 2 (ζ)[xe y -ye x ] , σ 2 = σ E (x, y)E E + σ T (x, y)E T . ( 23
)
As it will be shown in the next section, the lowest order macroscopic problem is a 1D beam problem, with extension and torsion. It thus involves macroscopic beam stresses which are simply defined from the integration over the crosssection S of the local or microscopic stresses σ 1 . Consequently the axial force T and the torque M take the form:

T (ζ) = S σ 2 ζζ dS , M (ζ) = S (-yσ 2 xζ + xσ 2 yζ )dS , (24)
and from the solution of problem ( 21), one can define the overall beam behavior such that:

   T M    = [k hom ]    E E E T    , (25)
where [k hom ] is the stiffness matrix, which is symmetric.

Macroscopic problem

The lowest order macroscopic problem can be derived from compatibility conditions which express that problem (21) admits a solution, see e.g. Buannic andCartraud (2000, 2001a). It amounts to integrate equilibrium equations of problem ( 21) over the cross-section. This yields:

                           dT /dζ = 0 dM /dζ = 0    T M    = [k hom ]    E E E T    E E = ∂u 1 ζ /∂ζ E T = ∂ϕ 1 /∂ζ , ( 26
)
with boundary conditions at ζ = 0 and ζ = εH. Since we have stress data for the 3D initial problem, and taking into account the overall equilibrium, these boundary conditions can be written as:

               T (0) = S t 0 • (-e Z )dS M (0) = S (yt 0 • e x -xt 0 • e y )dS T (εH) = T (0) M (εH) = M (0) , (27)
which corresponds to the application of the Saint-Venant principle, rigorously justified in the framework of asymptotic analysis of beams in Buannic and Cartraud (2001b).

The solution of this 1D macroscopic problem (26-27) is thus straightforward with a uniform macroscopic state: T = T (0), M = M (0), with the macroscopic strains E E and E T obtained from the inversion of (25) and u 1 ζ and ϕ 1 calculated thanks to (26) 4-5 and defined up to a constant.

Summary

One can summarize the results of the asymptotic expansion method with the following expressions:

u(x) = ε(u 1 ζ (ζ)e Z + ϕ 1 (ζ)[xe y -ye x ])+ ε 2 (χ E (x, y)E E + χ T (x, y)E T + u 2 ζ (ζ)e Z + ϕ 2 (ζ)[xe y -ye x ]) + O(ε 3 ) , σ = ε(σ E (x, y)E E + σ T (x, y)E T ) + O(ε 2 ) . (28)
It is recalled that microscopic fields χ E (x, y), χ T (x, y), σ E (x, y) σ T (x, y) are provided by the solution of the 2D microscopic problem (21) posed on the crosssection. Then, the expansions given in ( 28) can be easily computed up to the second-order rigid body motion, combining the previous solution of the 1D macroscopic problem with these microscopic fields.

Finite element solution

The variational formulation of the 2D microscopic problem (21) in the twisted coordinate system takes the form:

∀δu 2 (x, y), S ∇ s xy (δu 2 ) : σ 2 dxdy = 0, (29)
and from Eq. ( 21) 3 :

σ 2 = C : (∇ s xy (u 2 ) + ǫ macro ), (30)
with ǫ macro = ∇ s ζ (u 1 ). Hence one has:

∀δu 2 (x, y),
We recall that {∇ s xy (u 2 )} = L xy {u 2 }, see (19). Then a finite element approximation of the form {u 2 } = [N e ]{U e } is introduced , where [N e ] is the matrix of shape functions, and {U e } the nodal displacements, with three degrees of freedom at each node. The variational formulation yields:

[K]{U } = {F }, [K e ] = S e [N e ] T L T xy [C]L xy [N e ]dxdy, {F e } = - S e [N e ] T L T xy [C]{ǫ macro }dxdy, (32)
with [K] the stiffness matrix obtained from the assembly of element stiffness

matrices [K e ].
Note that in (32) the external load is given under the form of a macroscopic strain {ǫ macro }.

Once this system is solved, the stresses are computed thanks to (30) and after integration over the cross-section, the macroscopic beam stresses, i.e. the axial force and the torque are computed, thus providing the overall behavior [k hom ].

Validation of the homogenization approach

In this section, the microscopic response is computed for helical springs and seven-wire strands under axial loading. The 2D FE model based on helical homogenization has been implemented in an in-house code. This model is first validated for helical springs by comparison with an analytical solution. Another validation is also presented for seven-wire strands with a reference solution obtained from a 3D FE model.

For helical single wire or multi-wire structures subjected to a given macroscopic extension E E (E T = 0), first the 2D model is generated. The crosssection is meshed, with six-node triangle elements to improve the geometrical description as well as results accuracy. The solution of the microscopic 2D problem is defined up to a rigid body displacement in the twisted coordinate system, see Eq. ( 20), which can be fixed by prescribing the axial displacement u Z of an arbitrary node and the binormal displacement u y of a node on the line y = 0.

Then Eq. ( 32) is solved, and in the post-processing step, the computation of the axial force T and moment M are performed as well as the overall behavior.

Helical single wire structures

A helical single wire structure with circular cross-section is studied. R, Φ, n and a denote the helix radius, helix angle, number of helix pitches and the wire radius, respectively. Two types of structures can be distinguished: helical springs (large helix angle Φ and ratio R/a) and civil engineering cable (small angle Φ). The homogenization approach proposed in this paper is valid for any type of helical structures. However, in the literature, analytical solution is available only in the case of helical spring. Therefore, the validation of the homogenization approach is performed in that case.

The analytical solution may be found in [START_REF] Ancker | Pitch and curvature corrections for helical springs[END_REF]. When one end-section is clamped while the other is subjected to axial load T with a fixed rotation, the axial deflection δ at its end is given by:

δ = 4T R 3 n Ga 4 Ψ, Ψ = 1 - 3 16 ( a R ) 2 + 1 (1 + ν) tan 2 Φ ( 1 -ν 2 - ν 2 1 + 3 -7ν -20ν 2 -8ν 3 48(1 + ν) ( a R ) 2 + 1 + ν tan 2 Φ ) + ..., (33)
where Ψ is a pitch and curvature correction factor. The inputs of the analytical solution are the ratio a/R, the helix angle Φ and the Poisson coefficient ν. For given geometric and material parameters, Eq. ( 33) is used to compute the correction factor Ψ.

The numerical results provided by the homogenization approach are compared with the analytical solutions for helical springs as follows. For a given δ, the macroscopic strain E E = δ/nL, E T = 0 is imposed as the loading in (32) on the 2D FE model. Then the axial force T is computed. This leads to a numerical value of Ψ according to Eq. ( 33) 1 , which is compared to the analytical solution given by Eq. ( 33) 2 . For ν = 0.3, Fig. 4 shows the variation of the correction factor Ψ as a function of a/R for helix angle Φ = 70 • , 75 • , 80 • and 85 • . Only small differences between numerical and analytical results can be seen for a/R ≤ 0.2. This difference increases with a/R and as Φ decreases but remains less than 0.7% for Φ = 70 • and a/R = 0.35, which is small.

The same evolution of the differences between the numerical results and the analytical solution was observed in [START_REF] Jiang | Statically indeterminate contacts in axially loaded wire strand[END_REF], using a 3D FE model, with a free rotation. They are due to the non validity of the analytical model for large a/R and small helix angle Φ. However, our numerical results are in good agreement with those obtained from the analytical model providing a first validation of the computational homogenization approach.

Seven-wire strands

Multi-wire cables form a large class of civil engineering components. Sevenwire strands, composed of one layer of helical wires wrapped around a central wire, are the basic element of these cables. The major advantage of the twisted structure is its ability to carry large loads.

The static behavior of seven-wire strands was studied among others in Ghoreishi et al. ( 2007) using a 3D FE model. In that paper the overall strand stiffness was identified from computations performed on a model of two pitch length, and these results are considered as reference results in the following.

The static behavior is computed using the computational homogenization approach and the 2D FE model. The 2D mesh is generated as follows: an independent mesh for each wire of the seven-wire strand is first considered.

As mentioned before, the contact condition between the central and peripheral wires are assumed stick. Linear relations are then imposed at the contact point between the central and the peripheral wires, expressing the displacement continuity (u c = u p ), where the subscripts c and p correspond to the central and peripheral wires, respectively. In practice, the system (32) is condensed to take into account these conditions.

As an example, Fig. 6 shows the mesh of the cross-section of a strand with the following parameters: central wire with radius a and helical wires with helix radius R/a = 1.967 and angle Φ = 7.9 • . The cross-section of the central straight wire is circular. As for the previous helical single wire structure, the cross-section of helical wires is no longer circular in the (e X , e Y ) plane. Note that the helix radius R must be smaller than 2a, otherwise the adjacent helical wires would overlap each other. This example will be considered later in this section as well as in Part 2. Overall the macroscopic behavior of the seven-wire strand computed by the 2D FE model according to the homogenization approach is in good agreement with that obtained from the 3D model. This provides a second validation of the helical computational homogenization approach.

Lastly, microscopic displacements computed using the 2D FE model are analyzed. From the symmetry between the six helical wires, the displacements of only one peripheral wire is discussed. Fig. 9 shows the microscopic displacements in the cross-section of the seven-wire strand considered in Fig. 6 (R/a = 1.967, Φ = 7.9 • ), subjected to E E = 0.6%. The in plane component u 2

x of the central and the peripheral wire are shown in Fig. 9(a) and (c), respectively. One can observe the Poisson effect, with a linear evolution over the cross-section of u 2

x in the central wire, and an affine evolution in the peripheral wire, which is maintained in contact with the central wire. The axial displacement is presented in Fig. 9(b) and (d) for the central and the peripheral wire, respectively. One can notice that for the central wire the microscopic axial displacement is close to zero except in the vicinity of contact points where small variations occur. In the helical wire, a linear evolution of the microscopic axial displacement over the cross-section is found, due to local bending. For this example, the core wire radius is a = 2.7mm (the helical wire radius being 0.967a).

Material properties are: ν = 0.28 and E = 2.17e11P a. The computed axial force and torque are T = 190.3kN and M = 118.1N.m. This example will used in Part 2 of this paper, for wave propagation analysis in prestressed strands.

Conclusions

In this paper, the asymptotic expansion method has been applied to helical structures subjected to axial loads (traction and torsion) at its end sections.

Thanks to the use of a twisted coordinate system, the 3D elastic problem has been reduced to a 2D microscopic problem posed on the cross-section and a 1D macroscopic beam problem, which has an analytical solution. Therefore the main contribution of this work is the derivation of the 2D microscopic problem, which fully exploits the translational invariance of the problem. The solution of this problem enables the computation of the overall beam stiffness as well as microscopic stresses corresponding to a given macroscopic loading. The proposed approach has been validated for helical single wire structures and seven-wire strands and compares favorably with reference analytical results or 3D FE computations.

In Part 2 of this paper, the solution of the microscopic problem is used in order to take into account effects of prestress and geometry deformation on wave propagation.

Figure 1 :

 1 Figure 1: Left: One helix pitch and its twisted basis associated to the twisted coordinate system (x, y, Z). Right: view normal to the Z-axis. The point Z = s = 0 lies in the (e X , e Y ) plane.

Figure 2 :

 2 Figure 2: Cross-section of helical wires, R/a = 2 and (a) Φ = 10 • , (b) Φ = 30 • , (c) Φ = 50 • , (d) Φ = 70 • .

Figure 3 :

 3 Figure 3: 3D helical structures. (a) single wire, (b) seven-wire strand.

ζ

  = εZ denotes the slow scale or macroscopic 1D-variable and {x, y} denote the fast scale or microscopic 2D-variables. According to this change of variables, the differential operators become ∇ s (.) = ∇ s xy (.) + ε∇ s ζ (.) , ∇ • (.) = ∇ xy • (.) + ε∇ ζ • (.) , (16) where ∇ s ζ (.) and ∇ ζ • (.) correspond to partial differentiations with respect to the macroscopic variable ζ. ∇ s xy (.) and ∇ s xy •(.) denote the differential operators with respect to the microscopic variables x and y.

  2 ) : C : ǫ macro dxdy.

Figure 4 :

 4 Figure 4: Correction factor Ψ vs. a/R for Φ = 70 • , 75 • , 80 • , 85 • .

Figure 5 :

 5 Figure 5: Dimensionless microscopic displacements in the cross-section of a helical spring (R/a = 10, Φ = 75 • ) under axial deformation E E = 40%. (a) u 2 x /a, (b) u 2 ζ /a

Figure 6 :

 6 Figure 6: Mesh of seven-wire strand (R/a = 1.967, Φ = 7.9 • , a is the radius of the central wire)

Figure 7 :

 7 Figure 7: Dimensionless axial stiffness of seven-wire strand. k 11 vs. Φ. R/a = 2.

Figure 8 :

 8 Figure 8: Dimensionless stiffness coupling term of seven-wire strand. k 21 vs. Φ. R/a = 2.

Figure 9 :

 9 Figure 9: Dimensionless microscopic displacements of a seven-wire strand under axial deformation E E = 0.6%. (a) u 2 x /a and (b) u 2 ζ /a in the central wire (c) u 2 x /a and (d) u 2 ζ /a in the peripheral wire.


Segmentation of mesoscale ocean surface dynamics using satellite SST and SSH observations

Pierre Tandeo * , Bertrand Chapron † , Sileye Ba * , Emmanuelle Autret † , Ronan Fablet * * Institut Mines-Telecom, Telecom Bretagne, CNRS UMR 6285 LabSTICC, Pôle CID, Technopôle Brest Iroise, 29238 Cedex, France † IFREMER, Laboratoire d'océanographie spatiale, 29280 Plouzancé, France Abstract-Multi-satellite measurements of altimeter-derived Sea Surface Height (SSH) and Sea Surface Temperature (SST) provide a wealth of information about ocean circulation, especially mesoscale ocean dynamics which may involve strong spatiotemporal relationships between SSH and SST fields. Within an observation-driven framework, we investigate the extent to which mesoscale ocean dynamics may be decomposed into a mixture of dynamical modes, characterized by different local regressions between SSH and SST fields. Formally, we develop a novel latent class regression model to identify dynamical modes from joint SSH and SST observation series. Applied to the highly dynamical Agulhas region, we demonstrate and discuss the geophysical relevance of the proposed mixture model to achieve a spatiotemporal segmentation of the upper ocean dynamics.

Index Terms-SST, SSH, Observation-driven model, Latent class regression, Mesoscale ocean surface dynamics

I. INTRODUCTION

I N the last two decades, multi-satellite measurements of altimeter-derived Sea Surface Height (SSH) and multisensor measurements of Sea Surface Temperature (SST) have provided a wealth of information about ocean circulation and atmosphere-ocean interactions. As a depth-integrated quantity dependent upon the density structure of the water column, altimeter SSH estimations capture mesoscale structures, horizontal scales of 50 km to few hundred kilometers, and allow for the retrieval of surface currents using the geostrophy balance. This emerging and rich mesoscale circulation further stirs the large-scale SST fields. Accordingly, our picture of upper ocean dynamics has considerably evolved towards a complex system characterized by strong interactions, whose spatio-temporal variability extends over a wide range of scales. Furthermore, several studies (cf. [START_REF] Lacasce | Estimating subsurface horizontal and vertical velocities from sea-surface temperature[END_REF], [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF], [START_REF] Klein | Diagnosis of vertical velocities in the upper ocean from high resolution sea surface height[END_REF], [START_REF] Hausmann | The observed signature of mesoscale eddies in sea surface temperature and the associated heat transport[END_REF]) rationalize and demonstrate that fields of SST can become an active tracer coupled to the dynamics leading to strong correlations with SSH fields.

Such a framework can possibly guide the investigation and implementation of improved statistical means to optimally combine existing multi-altimeter SSH measurements with other satellite medium to high resolution observations (e.g., microwave sea surface temperature and salinity, scatterometer winds), augmented by the growing available in situ data (e.g., [START_REF] Ali | Estimation of ocean subsurface thermal structure from surface parameters: A neural network approach[END_REF], [START_REF] Chu | Determination of vertical thermal structure from sea surface temperature[END_REF], [START_REF] Wu | Estimation of subsurface Temperature anomaly in the North Atlantic using a self-organizing map neural network[END_REF]). Theoretically, the upper ocean turbulence for the horizontal scales between 50 km to few hundred kilometers is still consistent with the geostrophy turbulence theory. Under this assumption, the upper ocean dynamics may be simply predicted from surface density horizontal variations possibly dominated by SST variations. For such a case, a linear transfer function shall be identified between SSH and SST fields to also possibly lead to the estimation of the subsurface flow (e.g., [START_REF] Isern-Fontanet | Potential use of microwave sea surface temperatures for the estimation of ocean currents[END_REF], [START_REF] Lacasce | Estimating subsurface horizontal and vertical velocities from sea-surface temperature[END_REF]). This linear transfer function does not involve temporal differencing as in the maximum cross-correlation technique or alternate strategies (e.g., [START_REF] Bowen | Extracting multiyear surface currents from sequential thermal imagery using the maximum cross-correlation technique[END_REF], [START_REF] Marcello | Motion estimation techniques to automatically track oceanographic thermal structures in multisensor image sequences[END_REF], [START_REF] Chen | Surface velocity estimation from satellite imagery using displaced frame central difference equation[END_REF]). Note that other recent papers use nonlinear transfer functions to relate SST and SSH fields (e.g., [START_REF] Gaultier | On the inversion of submesoscale tracer fields to estimate the surface ocean circulation[END_REF]).

This strongly advocates for observation-driven studies to explore and characterize the local relationships between SST, SSH and the derived surface currents from satellite-based routine observations. Yet, as illustrated in Fig. 1, a simple linear transfer function cannot be expected to solely govern the whole mesoscale dynamics in a particular ocean region. As revealed, an overall spatial correlation exists, but for instance, relationships between SST gradients and altimetry-derived surface currents may spatially differ. In the warmer SST frontal zone, SST gradients correspond to large surface currents (top of the image). In the colder frontal area, large SST gradients do not reflect in large surface currents (bottom of the image). Besides, the clearly detected eddy (top-left of the image) is associated with weak SST gradients.

Within an observation-driven framework, one may consider joint PCA/EOF (Principal Component Analysis, Empirical Orthogonal Functions) procedures to decompose the relationships between SST and SSH fields, as widely used in ocean sensing applications (cf. [START_REF] Leuliette | Coupled pattern analysis of sea surface temperature and TOPEX/Poseidon sea surface height[END_REF], [START_REF] Casey | Sea surface temperature and sea surface height variability in the North Pacific Ocean from 1993 to 1999[END_REF]). Such EOF-based schemes would however resort to a single linear and global model. As such, this model could not address the spatial non-stationarity of the SST-SSH relationships illustrated in Fig. 1. By contrast, we here consider local linear transfer functions. We assume that locally, upper ocean dynamics may be analyzed according to a finite mixture model, where each component of the mixture is characterized by a local SST-SSH linear transfer function. This mixture-based representation relates to a nonlinear model with assumptions in accordance with the observations made above between SST and surface currents. In this paper, we propose to investigate such a model to (i) develop a probabilistic learningbased setting for the inference of such mixture models and the spatio-temporal segmentation of the identified dynamical modes (i.e., the different components of the mixture model), and (ii) evaluate the extent to which such mixture models are geophysically relevant to characterize the upper ocean dynamics over active ocean regions.

Hereafter we consider the Agulhas region, and the paper is organized as follows. Section II presents the remote sensing data and describes our probabilistic learning-based model. In Section III, the application to satellite observations is evaluated. We further discuss and summarize the key results of our investigations in Section IV.

II. DATA AND METHODOLOGY

A. Patch-based approach

As mentioned above, recent theoretical and numerical experiments have stressed that upper ocean dynamics may be characterized by couplings between SSH and SST according to the following relationship in the Fourier domain (cf. [START_REF] Held | Surface quasi-geostrophic dynamics[END_REF]):

F H ( SSH) = -γ|k| -α F T ( SST) ( 1
)
where k is the horizontal wavenumber vector, F H and F T are linear filters of SSH and SST respectively. The α parameter sets up the effective coupling between surface fields. For α = 1, Eq. ( 1) resorts to the Surface Quasi-Geostrophy (SQG) model. In [START_REF] Isern-Fontanet | Potential use of microwave sea surface temperatures for the estimation of ocean currents[END_REF], F H and F T were band-pass filters between 80 km and 300 km. As α increases, the smoothing increases and couplings decreases. For α = 2, the filtered SST would trace the vorticity. Formally, Eq. ( 1) states that surface currents can be regarded as spatial derivatives of a filtered version of the SST field. The parameter γ relates to a normalization constraint. In general, parameters γ and α as well as the definition of the filters F H and F T , may spatially vary such that a single linear transfer function as in Eq. ( 1) is unlikely to apply globally as illustrated in Fig. 1. These considerations led us to hypothesize that zonal and meridional surface currents noted (U,V) and SSH can still locally relate to SST derivatives, but according to a finite set of K linear transfer functions, hereafter referred to as K dynamical modes. Formally, this is stated in the Fourier domain as

( SSH, U, V) = H k ( SST) (2)
where H k characterizes the k th dynamical mode, which locally relates SSH and SST fields through linear filter H k . In this study, we do not consider any band-pass filters F H and F T . Using a matrix formulation, Eq. ( 2) is rewritten in the real domain as a patch-based linear regression

Y(s i , t i ) = H k (X(s i , t i )) (3)
where Y(s i , t i ) encodes the local SSH variability through a 3-dimensional vector formed by the SSH value and the surface current (U,V) at spatio-temporal location (s i , t i ) and X(s i , t i ) is the vectorized version of the local SST patch centered in s i at time t i (cf. Fig. 2). It may be noted that we encoded local SSH variations at spatio-temporal location (s i , t i ) through the surface currents which are computed as the spatial derivatives of the SSH field. As such, it constrains the method to account for spatial regularity. The linear operator associated with dynamical mode k is corresponding to the local version of H k in Eq. ( 2). It corresponds to three vectorized versions of spatial convolution matrices. Here, p defines the size of the local SST neighborhood around s i and is set according to the Rossby radius of the study region, i.e. the mean size of the mesoscale ocean structures like eddies. For the Great Agulhas current region, we set it up to 200 km, i.e. p = 81 for the spatial resolution of the considered data.

B. Remote sensing data

As SSH and surface geostrophic current (U,V) data, we use the daily delayed time Maps of Absolute Dynamic Topography (MADT) produced by Collecte Localisation Satellites (CLS) available online at http://www.aviso.oceanobs.com/. This information combines the signal of several altimeters onto a 1/3 • Mercator projection grid. We use the 2004 data since four altimeters were available (Jason-1, Envisat or ERS-2, Topex/Poseidon and GFO). As SST data, we use optimally interpolated microwave SSTs provided by Remote Sensing System (RSS) available online at http://www.ssmi.com/. It combines the signal of three microwave radiometers (TMI, AMSR-E and WindSAT) which are robust to the presence of clouds. The spatial resolution is 1/4 • × 1/4 • and the temporal

X(s i , t i ) Y (s i , t i ) √ p X(s j , t i ) Y (s j , t i ) X(s i , t i ) Y (s i , t i ) X(s j , t i ) Y (s j , t i ) Fig. 2.
Sketch of SST patches (in degree represented in false colours, the more red the hoter), noted X, and the corresponding SSHs (in meter represented by dots, the bigger the higher) and surface currents (in meter per second represented by quivers, the longer the stronger) noted Y at the central location s i and s j at time t i . resolution is the same as the MADT data, i.e. daily. We bilinearly interpolate the MADT data onto the SST grid. We focus on the Agulhas region between longitudes 5 • E to 65 • E and latitudes 30 • S to 48 degreeS.

Given the joint series of satellite observations, we extract SST patches (noted X) and the associated SSH with surface current (U,V) at the center of the patches (noted Y). Overall, the processed dataset is composed of ∼ 5×10 6 pairs of vectors X and Y. To infer the parameters of the considered mixture model, i.e. the parameters of each dynamical mode in Eq.

(3), we first build a training dataset as a random sample of n = 10 5 elements (for a given day, we use about 2% of the data to fit the model). In a second step, we apply the inferred mixture model to the entire processed dataset to extract the spatio-temporal of the different dynamical modes.

C. Latent class regression model

Our objective is to identify K different dynamical modes (latent variable Z) from a joint set of SST patches (pdimensional vector X) and SSH with zonal and meridional surface currents (3-dimensional vector Y). In this paper, we assume that the conditional probability of Y given X and the dynamical mode Z = k is given by

p (Y|X, Z = k) ∝ N k (Y; Xβ k , Σ k ) (4)
where N k represents a multivariate Gaussian probability density function evaluated in Y with mean Xβ k and covariance Σ k . Hence, the conditional probability of Y|X resorts to a mixture of Normal distributions

p (Y|X, θ) = K k=1 λ k N k (Y; Xβ k , Σ k ) (5)
where λ k is the prior probability of mode k. To simplify the notations, we store the overall parameters of Eq. ( 5) in

θ = (λ 1 , . . . , λ K , β 1 , . . . , β K , Σ 1 , . . . , Σ K ).
In the literature, this model is referred to as a "latent class regression" or "clusterwise regression" (cf. [START_REF] Desarbo | A maximum likelihood methodology for clusterwise linear regression[END_REF]). By construction, it imposes that 0 ≤ λ k ≤ 1, K k=1 λ k = 1 and Σ k is positive defined. The maximum likelihood estimation procedure for model parameters θ is given below.

D. Model learning

To learn model parameters θ in Eq. ( 5), we resort to a maximum likelihood criterion and use an iterative Expectation-Maximization (EM) procedure (cf. [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]). It relies on the maximization of the log-likelihood given by

L (θ) = n i=1 log (p (Y(s i , t i )|X(s i , t i ), θ)) ( 6
)
where n is the number of observations of the training dataset described in Section II-B. From a given initialization, the EM procedure iterates an E-step (Expectation-step) and M-step (Maximization-step). At a given iteration, using the Bayes theorem, the E-step resorts to the computation of the posterior probabilities of the latent variable Z for each spatio-temporal location (s i , t i ) given current parameter estimate θ:

πk (s i , t i ) = λk N k Y(s i , t i ); X(s i , t i ) βk , Σk p Y(s i , t i )|X(s i , t i ), θ , ∀k. (7)
The M-step then minimizes the expectation of the loglikelihood conditionally to the current parameter estimate θ. This leads to the update of the prior probabilities of the latent variable Z as

λk = n i=1 πk (s i , t i ) n , ∀k. (8)
The updated regression parameters βk , ∀k are derived by fitting K distinct linear regressions using a weighted least square criterion on the n observations where the weights are given by the posterior probabilities given in Eq. ( 7) as in [START_REF] Schlittgen | A weighted least-squares approach to clusterwise regression[END_REF].

Then, we maximize L with respect to Σ k and obtain

Σk = n i=1 πk (s i , t i )ǫ k (s i , t i ) ⊤ ǫ k (s i , t i ) n i=1 πk (s i , t i ) , ∀k (9)
where ǫ k (s i , t i ) = Y(s i , t i ) -X(s i , t i ) βk and ⊤ represents the transpose. The algorithm iterates the E-step and M-step until a negligible increase of the log-likelihood L which is strictly growing.

A key aspect of the latent class regression is the choice of K, the number of clusters. Different statistical criteria state the selection of parameter K as a trade-off between the likelihood and the complexity of the model (cf. [START_REF] Kadane | Methods and criteria for model selection[END_REF]) such as the Akaike Information Criterion (AIC). However, the optimization of these criteria makes no effort to distinguish the error explained by the regression fit and the error explained by the clustering process. In practice, there is an actual potential for overfitting with latent class regression model (cf. [START_REF] Brusco | Cautionary remarks on the use of clusterwise regression[END_REF]). Then, for a given number of clusters K, we suggest different evaluations of the EM algorithm to reach a greater consistency in the estimation of model parameters (cf. [START_REF] Hennig | Identifiablity of models for clusterwise linear regression[END_REF]). The idea is to use random values πk as initialization values of the EM procedure and select parameter estimates corresponding to the greatest likelihood (see [START_REF] Biernacki | Choosing starting values for the EM algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF] for more details). In our case, the inference based on AIC would lead to K values between 4 and 9. We performed a complementary qualitative analysis and the setting K = 4 resulted as a good trade-off between the geophysical interpretation of the model and regression error statistics (i.e., the maximization of the likelihood).

E. Spatio-temporal segmentation of dynamical modes and SSH/current predictions

We exploit the inferred mixture model with parameter θ to perform a spatio-temporal segmentation of the underlying dynamical modes. More precisely, for any spatial location s and time t within the Agulhas region over the year 2004 (cf. processed data in Section II-B), we use Eq. ( 7) to evaluate the posterior probability πk (s, t) for the K = 4 dynamical modes. Then, the pixels will be assigned to the most likely dynamical mode. One can also estimate for each time t, the relative spatial occurrence of each dynamical mode using Eq. ( 8). Besides, using Eq. ( 3), the estimation of the SSH and surface current at the spatial location s and time t from the associated SST patch follows from the fuzzy regression

Ŷ(s, t) = K k=1 πk (s, t)X(s, t) βk (10)
where the βk parameterize the linear operator H k . They correspond to the p×3 regression coefficient matrices between the patch of SSTs (X in degree) and the central value of SSH along with zonal U and meridional V components of the surface currents (Y in meter and meter per second). All the locations of the Agulhas current share the same K matrices βk which are not changing temporally nor spatially. It may be outline that no additional constraint is set on the posterior probabilities πk (s, t), which could reveal space-time variations of the distribution of the dynamical modes (including seasonal variations).

III. CHARACTERIZATION OF OCEAN SURFACE DYNAMICS

We first report the temporal evolution of the relative spatial occurrence of the K = 4 dynamical modes (cf. Fig. 3). The dynamical modes involve clear seasonal cycles. Dynamical modes 1 (red) and 4 (blue) depict similar temporal variations, completely out of phase with dynamical modes 2 (green) and 3 (cyan). To study the spatial distribution of these dynamical modes, we focus on two dates corresponding to the maximal and minimal values of the seasonal cycle, namely the 1 st of March and the 1 st of September, 2004. They correspond to the maximum or minimum values of the λk .

For these two dates, from the maps of posterior probabilities πk (s, t), we determine the segmentation maps of each dynamical mode as illustrated in Fig. 4. The animations of the time series of these daily maps in the Agulhas current over 2004 are available as supplementary material or online at: http://tandeo.wordpress.com/communications/articles/. A first qualitative analysis of these maps highlights a clear spatiotemporal segmentation of the different dynamical modes that Fig. 3. Time series of the relative proportion of locations associated with the four dynamical modes within the Agulhas region for the year 2004. In the subsequent, colors red, green, cyan and blue respectively distinguish the first, second, third and fourth dynamical modes. The straight lines correspond to a six month time separation: 1 st of March and 1 st of September, 2004. The corresponding spatial maps of the posterior probabilities for each dynamical mode are given in Fig. 4.

can be interpreted from a geophysical perspective in terms of different geophysical processes. We also report for each dynamical mode the observed distributions of current, height and temperature values (cf. Fig. 5). The first dynamical mode (red) characterizes very strong current magnitude and warm waters. It is primarily associated with the main Agulhas current that flows down the East coast of Africa through the Agulhas ridge. This mode also involves mesoscale eddies, the so-called warm core Agulhas rings with strong surface currents, low temperature gradients and middle-range SSH values around 0.5 m. The later seems to be a discriminative feature of this first mode. The second dynamical mode (green) mainly relates to the eastward Agulhas return current that hits a part of the South Atlantic current. It creates a subtropical front varying from 36 • S to 44 • S with strong eastward currents, middlerange SST gradients and large SSH values (about 1 m). This front was clearly observed in the upper part of Fig. 1(b). The third (cyan) and fourth (blue) dynamical modes correspond to weaker surface currents. The third one is characterized by mid-temperatures and westward currents whereas the fourth one involves colder temperatures and eastward currents. Let us stress that the third dynamical mode involves large SST gradients but weak surface currents as identified in the lower part of Fig. 1(b). In this mode, the SST can be clearly identified as a passive tracer of the surface upper ocean dynamics.

To characterize more precisely the inferred model, we plot the regression line and the 95% confidence region from the estimated parameters βk and Σk of each dynamical mode (cf. Fig. 6). These results clearly stress the relevance of a mixture model, compared to a single linear model (represented as the black line). Similar slopes are observed for the second and third (green and cyan) dynamical modes as well as the first and fourth (red and blue) ones. Stronger differences among the dynamical modes are outlined regarding the associated regression error variances. We can notice that the the first mode involves the greater variance and the fourth one the lower one.

We further investigate the geophysical consistency of the identified dynamical modes from the comparison between the true MADT SSH and surface currents in Fig. 7(a) and the prediction of the latent class regression model in Fig. 7(b) given by Eq. [START_REF] Desarbo | A maximum likelihood methodology for clusterwise linear regression[END_REF]. Overall, as reported in Table I, a good agreement is obtained with a global correlation coefficient of 0.72 for the surface currents and 0.96 for the SSH; it can locally be very large as illustrated in the right column of Fig. 7 (corresponding to the zone depicted in Fig. 1). This zone involves the four dynamical modes. The mixture model enables us to retrieve both the large warm eddy associated with weak SST gradients (top-left of the zone and first dynamical mode), the relatively large surface currents along the large warmer SST gradients (upper part and second mode), as well as the rather weak currents along the large but colder SST front (lower part and third mode). For comparison purposes, we also plot the results issued from a single linear transfer function in Fig. 7(c). This model clearly underestimates the surface Fig. 6. SSH as a function of SST. For each dynamical mode, we give the regression line and the 95% confidence envelope. The four regressions are highly significant (p-values < 10 -3 ). The black line is a benchmark curve corresponding to the global linear regression (with a one-mode model).

currents within the warm eddy (top-left) and overestimates the currents of the colder frontal zone (lower part), which stresses the requirement for considering a mixture model.

To further characterize each dynamical mode, we report in Table II the correlation and RMSE statistics computed within the associated spatio-temporal domain, i.e. the domain comprising all spatio-temporal locations assigned to dynamical mode k according to posterior probabilities πk (s, t) computed in Eq. [START_REF] Chen | Surface velocity estimation from satellite imagery using displaced frame central difference equation[END_REF]. We compare the latent class regression model with respect to both the true MADT data and a SQG-like hypothesis, i.e. Eq. ( 1) with α = 1, within the space-time region associated with each dynamical mode. This analysis clearly discriminates the second and fourth modes from the first and third modes. The linear transfer functions of the second and fourth dynamical modes capture a large part of the variability of the true SSH data. These first two modes have also a good consistency with the SQG hypothesis, correlation coefficients of 0.63 and 0.68, respectively. These results suggest that the SST might be regarded as an active tracer of the surface dynamics in the associated regions. By contrast, the SQG hypothesis poorly fits to the first and third dynamical modes, with correlation coefficients of 0.33 and 0.27, respectively. These two dynamical modes also involve a slightly lower predictability of the linear transfer functions to retrieve the SSH and surface currents.

To summarize, the four dynamical modes correspond to different physical parameter values in Eq. ( 1). The second and fourth modes seem to correspond to α ≃ 1 (close to the SQG model) whereas the first and fourth modes appear to be α > 1 (coupling of SST and SSH at large scales). The factor γ relates to the amplitude of the surface geostrophic currents: large values relate to strong currents as retrieved for modes 1 and 2, and low values relate to strong currents as observed for modes 3 and 4. Overall, these results are consistent with the previous work reported in [START_REF] Isern-Fontanet | Potential use of microwave sea surface temperatures for the estimation of ocean currents[END_REF] and [START_REF] Xu | The effects of altimeter instrument noise on the estimation of the wavenumber spectrum of the sea surface height[END_REF]. In particular, Isern-Fontanet et al. explored the SQG hypothesis from a phase-correlation analysis between SSH and SST fields, while Xu et al. explored the coupling assumption from a spectral analysis of SSH fields. These authors concluded that SQG-like dynamics would mainly occur near the edge of the large current system. Compared to these analysis, our contribution is two-fold: the quantitative characterization of the extent to which the SQG dynamics applies through correlation statistics as well as the actual space-time tracking of the regions associated with SQGlike and non-SQG dynamical modes.

IV. CONCLUSION AND PERSPECTIVES

In this paper, we propose an observation-driven framework to identify, characterize and track ocean surface dynamical modes. We rely on a latent class regression model, where the dynamical modes are characterized by a local linear transfer function between SST, SSH and surface current (U,V), in agreement with the theoretical assumption given in Eq. ( 1). This probabilistic approach locally relates the distribution of the SSH and sea surface currents conditionally to the SST via a nonlinear model: a Gaussian mixture of linear transfer functions. The statistical parameters of the model are estimated using a maximum likelihood approach.

We applied the proposed methodology to the 2004 daily 1/4 • × 1/4 • satellite SST and SSH image series. The reported results retrieved a relevant spatio-temporal decomposition of ocean surface dynamics in the Agulhas region according to four dynamical modes: (i) the main Agulhas current and warm core rings characterized by strong currents and hot temperatures where the SST is weakly correlated with SSH, (ii) the return Agulhas current with lower temperatures and currents where the SST is an active tracer, (iii) local front regions where strong SST gradients do not seem to affect the current velocities, and (iv) a weaker dynamical mode where SST is strongly correlated to SSH.

Our study complements previous theoretical studies which showed that mesoscale upper ocean dynamics may be characterized by a linear coupling between SST and SSH (cf. [START_REF] Lacasce | Estimating subsurface horizontal and vertical velocities from sea-surface temperature[END_REF], [START_REF] Lapeyre | Dynamics of the upper oceanic layers in terms of surface quasigeostrophy theory[END_REF], [START_REF] Isern-Fontanet | Three-dimensional reconstruction of oceanic mesoscale currents from surface information[END_REF], [START_REF] Klein | Diagnosis of vertical velocities in the upper ocean from high resolution sea surface height[END_REF], [START_REF] Hausmann | The observed signature of mesoscale eddies in sea surface temperature and the associated heat transport[END_REF]). Following a fully observation-driven framework, the proposed latent regression model enabled us to identify different dynamical modes, including some SQGlike ones, and to track the space-time extension of each dynamical mode. The reported results clearly pointed out the requirement for considering a mixture model to decompose the space-time variabilities of the ocean surface dynamics. Regarding methodological aspects, it may be pointed out that EOF-based schemes (cf. [START_REF] Leuliette | Coupled pattern analysis of sea surface temperature and TOPEX/Poseidon sea surface height[END_REF], [START_REF] Casey | Sea surface temperature and sea surface height variability in the North Pacific Ocean from 1993 to 1999[END_REF]) could not reveal such nonstationary space-time variabilities. Joint EOF scheme typically decomposes a global linear mapping between the analyzed fields according to principal modes. In our case, such an EOF decomposition could be considered for each dynamical mode to further characterize the associated linear transfer function with respect to joint SST-SSH principal modes.

Regarding our future work, we will further investigate latent class regression models with additional regressors. Among others, it seems appealing to explore how time-lagged SST features and other geophysical fields such as wind speed, mixed-layer depth, salinity and chlorophyll-a concentration (cf. [START_REF] Saraceno | On the relationship between satellite-retrieved surface temperature fronts and chlorophyll a in the western South Atlantic[END_REF]) could improve the estimation of SSH and surface currents. We also plan to apply the proposed model to other strongly active ocean regions such as the Gulf Stream system. Our objective will be to determine shared and/or systemspecific dynamical modes. Future work will also investigate more detailed physical interpretation of the identified dynamical modes, especially in terms of spectral characteristics. For instance, it would be interesting to relate more precisely the physical parameters γ and α of Eq. ( 1) to the different hidden dynamical modes extracted by our statistical approach. Whereas factor γ seems to be well estimated in this present paper, the spatial resolution of the satellite data (up to 50 km) does not permit to detect dynamical modes with α < 1, corresponding to more local couplings between geophysical fields. Besides, such an improved model shall then possibly address both (i) the higher resolution prediction of mesoscale ocean surface currents from SST spatio-temporal fields (and improve the non geostrophic components estimation as in [START_REF] Marcello | Motion estimation techniques to automatically track oceanographic thermal structures in multisensor image sequences[END_REF], [START_REF] Chen | Surface velocity estimation from satellite imagery using displaced frame central difference equation[END_REF]), and (ii) the extraction of new local and global descriptors of ocean surface dynamics from satellite sea surface observations (cf. [START_REF] Ba | Statistical descriptors of ocean regimes from the geometric regularity of SST observations[END_REF], [START_REF] Fablet | Multiscale geometric deformations along planar curves: application to satellite tracking and ocean observation data[END_REF]).

  32

Fig. 5 .

 5 Fig. 5. Distribution of the SST (a), SSH (b), surface current (U,V) norm (c) and direction (d). The results are given for each dynamical mode within the Agulhas current for the whole year 2004.

(a) 1 Fig. 4 .

 14 Fig.4. Maps of the posterior probabilities given in Eq. (7) of the dynamical modes given the SST and SSH fields, the 1 st of March (a) and the 1 st of September (b), 2004 within the Agulhas current. We use a four-class latent regression model fitted from the whole year 2004 (see text for details). For a given location and time, the sum of the four probabilities is equal to 1. The animations of the daily maps is available as supplementary material.

Fig. 7 .

 7 Fig. 7. True Maps of Absolute Dynamic Topography (MADT) data (a) and the results for the proposed latent class regression model (b) using Eq. (10) and for the linear model (c), the 1 st of January, 2004 within the Agulhas current. The left column corresponds to the SSH and the right column corresponds to the SSH and surface currents (U,V) for the zone depicted in Fig.1and in the black box.

TABLE I CORRELATION

 I AND RMSE STATISTICS WITHIN THE AGULHAS CURRENT FOR THE WHOLE YEAR 2004. THE LABELS "MADT", "LINEAR" AND "LATENT" REFER RESPECTIVELY TO THE TRUE MADT DATA, THE LINEAR MODEL AND THE PROPOSED LATENT CLASS REGRESSION MODEL. ALL THE CORRELATIONS ARE STATISTICALLY SIGNIFICANT (P-VALUES < 10 -3 ).

	Correlation (RMSE)	SSH	(U,V)
	MADT LINEAR	0.89 (0.22) 0.61 (0.19)
	MADT LATENT	0.96 (0.16) 0.72 (0.16)

TABLE II CORRELATION

 II AND RMSE STATISTICS FOR SPATIO-TEMPORAL LOCATIONS ASSIGNED TO EACH DYNAMICAL MODE ACCORDING TO POSTERIOR PROBABILITIES COMPUTED IN EQ. (7) WITHIN THE AGULHAS CURRENT FOR THE WHOLE YEAR 2004. THE LABELS "MADT", "SQG" AND "LATENT" REFER RESPECTIVELY TO THE TRUE MADT DATA, A SURFACE QUASI-GEOSTROPHIC HYPOTHESIS AND THE PROPOSED LATENT CLASS REGRESSION MODEL. ALL THE CORRELATIONS ARE STATISTICALLY SIGNIFICANT (P-VALUES < 10 -3 ).

		1.5		36 o S
	3 2 o S	1		38 o S
	3 6 o S	0.5	SSH (m)	40 o S
	o S 4 0	0		42 o S
	o S 4 4 4 4 o S 4 0 o S 3 6 o S 3 2 o S 4 4 o S 4 0 o S 3 6 o S 3 2 o S	Correlation (RMSE) SSH MADT SSH LATENT (U,V) MADT (U,V) LATENT 24 o E 36 o E 48 o E 24 o E 36 o E 48 o E (b) Latent class regression Mode 1 Mode 2 0.72 (0.26) 0.93 (0.16) 0.97 (0.13) 0.96 (0.12) Mode 3 Mode 4 0.62 (0.33) 0.88 (0.15) 0.63 (0.08) 0.88 (0.06) 6 0 o E -0.5 32 o E 34 o E (a) True MADT data 6 0 o E SSH (m) -0.5 0 0.5 1 1.5 32 o E 34 o E 42 o S 40 o S 38 o S 36 o S 24 o E 36 o E 48 o E 6 0 o E SSH (m) -0.5 0 0.5 1 1.5 32 o E 34 o E 42 o S 40 o S 38 o S 36 o S (U,V) 12 o E 12 o E 12 o E (c) Linear regression	36 o E 36 o E 36 o E	38 o E 38 o E 38 o E	40 o E 40 o E 40 o E	42 o E 42 o E 42 o E

MADT (U,V) SQG 0.33 (0.41) 0.63 (0.28) 0.27 (0.25) 0.68 (0.18)

V. ACKNOWLEDGMENTS

We would like to thank the Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO) and the RSS projects for respectively providing the altimeter-derived SSH and surface current, and the microwave SST data. This study was partially supported by grant CREATE FI2O, the LabexMER, ANR RedHots and Institut Mines-Telecom. We are grateful to the reviewers and the editor for their numerous constructive comments.


Thiago Peixoto De Araújo
email: thiagopeixotoaraujo@gmail.com

Fernando Machado De Mello

Essaid Bilal
email: bilal@emse.fr

Resumo

Digital Image Processing Techniques for Landslides Characterization in Rio de

Keywords: Remote Sensing, IHS-RGB Transformation, Landslides, Rio de Janeiro, Environmental Geology

The use of images by orbital remote sensing has shown great utility in environmental geology and engineering by mapping surveys and researches. The advantage of employing such remote sensing technique is due not only to its wider coverage area but also to its high spatial resolution quality, allowing a range of possibilities by means of operations among electromagnetic bands. The main purpose of this paper is to do comparisons among Landsat images of TM and ETM+ sensors, in Rio de Janeiro Region, Brazil, by means of color compositions that come from IHS-RGB Transformation and compositions that come from Ratios Band as well. Our studies point that ETM+ images presented better results in landslides studies due to their better accuracy in comparison to TM images.

INTRODUCTION

In recent years there has been a growing demand for environmental geotechnological information applied to geology, geotechnics and other related areas, especially where it is intended to carry out projects for the implementation of conservation units, intervention or environmental restoration of degraded areas. The production of reliable GIS maps is an indispensable stage either to the supply of technical analysis or to planners to support important regional decisions, where they propose, for example, assessments of environmental sustainability of an enterprise that must take into account physical, economic and social, beyond environmental impact.

By integrating both alphanumeric and geographical data through surveys and satellite images, and even aerial photographs, it is possible to make corrections in those images, providing they have a radial distortion, because the center of an image perspective is not met in the infinity, thus giving a conical projection in response. By correcting these images, we can make several characterization studies of the environmental areas and even monitor these areas as well.

Remote sensing data by orbital sensors can be used for motion studies of gravitational mass movement in both direct detections of the consequences of these processes (scars, exposed rocks, landslides) and through the identification of indicators, such as: use of land, morphology and dissection of relief [START_REF] King | Spatial assessment of erosion: contribution of remote sensing, a review[END_REF][START_REF] Berger | Geoindicator checklist[END_REF]. Changes in spectral characteristics of surface, detected by sensorial systems, are the result of the vegetation removal and surface layers of soil by either natural or anthropogenic processes.

Differences in spectral behavior among soil, vegetation (type and size) and anthropogenic features are observed in the bands of visible and infrared data obtained by Landsat TM and ETM+ sensors. Differences between scars with soil or exposed rocks and scars with resumption of vegetation (grasses, ferns) can be detected by the sensors.

For Drury (1993) and [START_REF] Walsh | An overview of scale, pattern, processes relationships in Geomorphology: a remote sensing and GIS perspective[END_REF] processing such as Ratio Band, for example, highlights those differences, serving as indicators of the geomorphologic processes and geotechnical problems. [START_REF] Mckean | Remote sensing and landslide hazard assessment[END_REF] and [START_REF] Gao | Micro-scale modelling of terrain susceptibility to landsliding from a DEM: a GIS approach[END_REF] say that for landslides, the vegetation covers removal and consequent exposure of soil or rock and deposited transported material along the slope, producing variations of shade, color, tint and texture in image. Specific forms are observed, thus allowing direct interpretation of these features.

The shapes of the scars caused by landslides, according to [START_REF] Hansen | Strategies for classification of landslides[END_REF] and [START_REF] Patton | Drainage basin morphometry and floods[END_REF], are elliptical or conical and when new, are observed from the head to the base by sliding. These features usually develop either in concave slopes along the drainage line or close to sites with anthropogenic changes, such as roads, extending over high or medium slope and depositing the removed material in low valleys, where formation of talus occurs, altering the topography of the place [START_REF] Patton | Drainage basin morphometry and floods[END_REF][START_REF] Gao | Micro-scale modelling of terrain susceptibility to landsliding from a DEM: a GIS approach[END_REF][START_REF] Sestini | Variáveis geomorfológicas no estudo de deslizamentos em Caraguatatuba -SP utilizando imagens TM-LANDSAT e SIG[END_REF]. However, sites susceptible to landslides can be identified through the interpretation of direct and indirect elements in the orbital images.

So, the study and development of techniques of remote sensing through the operations IHS-RGB Transformation and Ratio Bands is justified because these techniques allow us to individualize targets/objects to be studied, by using electromagnetic band [START_REF] Araujo | Processamento de imagens digitais -Razão entre Bandas[END_REF].

However, the main objective of this study is to demonstrate not only the advantages but also the disadvantages of the techniques of IHS-RGB Transformation and Ratio Bands held on Landsat images. Through images and some indications, it is observed that they are more efficient for certain studies, such as landslides [START_REF] Sestini | Variáveis geomorfológicas no estudo de deslizamentos em Caraguatatuba -SP utilizando imagens TM-LANDSAT e SIG[END_REF], characterization of geoenvironmental areas, erosion and many others.

The area chosen for this study is in the State of Rio de Janeiro, southeastern Brazil, bordering the states of São Paulo, Minas Gerais and Espirito Santo. The State of Rio de Janeiro is famous all over the world for its peculiar natural beauty. However it has greatly suffered by environmental problems (landslides, for example) and many social problems as well. FIGURE 1 shows the location of Rio de Janeiro in the context of southeastern Brazil.

MATERIAL AND METHODS

In this work, satellite images by Landsat 5 with the sensor TM (Thematic Mapper) and Landsat 7 with the sensor ETM+ (Enhanced Thematic Mapper Plus) have been used. The scene of Rio de Janeiro state used in this study was 217/76. These images are available free of charge on the website of National Institute for Space Research -Brazil (INPE). According to the issued report and to the attached images, the features in the acquisition of the scenes are: Scene 217/76 -Datum SAD 69; SAD 69 reference ellipsoid; unit in meters; UTM cartographic system; zone -23 South; and Product Framing Method Path 217 and Row 076.

The methods used in this study were the operations of IHS-RGB Transformation and the Ratio Bands, both procedures performed in the program ENVI™ (Environmental for Visualizing Images).

IHS-RGB TRANSFORMATION

To make the operation of IHS-RGB Transformation, it is necessary to have the individualized electromagnetic bands and then scroll through the sequence of commands in the program main menu: "File -Open External File -Landsat -GeoTIFF". Afterwards load these bands. After that, select which band will be in R, G or B and load the image (in "Load RGB"). The transformation will take place and you will get a color composition as result.

BAND RATIOS

For Ratio Bands operation, the procedure for charging the electromagnetic bands is the same as previous, but the sequence of commands to be covered is as follows: "Transform -Band Ratios". Then we should indicate which band will be in the numerator and which band will be on the denominator, then click "Enter Pair" and "ok". After that, you must provide a name and place of exit and "ok". Then you should proceed as done in IHS-RGB Transformation, in other words, indicate what reason will be in R, G or B and click "Load RGB".

RESULTS

When processing a satellite image in color composite, a transformation IHS-RGB is held. IHS are the initials of intensity (I), hue (H) and saturation (S), where I, or brightness, is the measure of total energy involved in all the wavelengths, being responsible for the sensation of brightness of the incident energy on the eyes; H, or color of an object is a measure of the average wavelength of light that is reflected or issued, defining the object color; and S, or purity, expressed the range of wavelength around the mean wavelength at which energy is reflected or transmitted. A spectrally pure color is the result of a high value of saturation, whereas low value of saturation indicates a mixture of wavelengths producing pastel tones.

However, the transformation is done through the mathematical algorithms, performed by the program, which lists the RGB space to IHS space. The equations 1, 2 and 3 show these mathematical algorithms used to make the transformation.

 

B G R I    3 1 (1)                                   2 1 2 1 5 , 0 cos B G G R G R B R G R H (2)       B G R B G R S , , min 3 1     (3)
The electromagnetic bands are nothing else but a portion of the electromagnetic spectrum bounded by two wavelengths.

Landsat 5 shows the scanner TM (Thematic Mapper), with features such as 7 bands with a spatial resolution of 30 meters, in the infrared channels (channels 1, 2, 3, 4, 5 and 7), and 120 meters, in thermal channel (band 6). As for Landsat 7 the channels of the visible and infrared regions (bands 1, 2, 3, 4, 5 and 7), remained with the same spatial resolution of 30 meters. However, the thermal infrared bands, now in high (H) and low (L) frequency (bands 6H and 6L), have a resolution of 60 meters, against 120 meters, the thermal channel (band 6), in Landsat 5.

And a major advance in the ETM + of Landsat 7, was the channel panchromatic (band 8), offering a spatial resolution of 15 meters. On TABLE 1 it is illustrated a comparative table of maximum and minimum limits of wavelengths (μm) of spectral bands of TM and ETM + sensors by their respective satellites, Landsat 5 and Landsat 7. A mathematical operation Ratio Bands discriminates subtle differences in spectral behavior of different targets because in the original bands only gross differences are observed. Thus, observing the spectral behavior of targets of interest for the application of ratio bands is selected aiming at their maximum and minimum values of reflectance, so that they express the gradients of the spectral curve of the objects of interest giving the highlight of these targets [START_REF] Sestini | Variáveis geomorfológicas no estudo de deslizamentos em Caraguatatuba -SP utilizando imagens TM-LANDSAT e SIG[END_REF].

According to some authors, the Ratio Bands has advantages [START_REF] Harrison | Introduction to remotely sensed data[END_REF]and Envi Guide, 2006) and disadvantages [START_REF] Mather | Computer processing of remotely sensed images: an introduction[END_REF][START_REF] Crosta | Processamento digital de imagens de sensoriamento remoto[END_REF]Drury, 1993).

The main advantages are: (i) The difference in spectral response of the same target in different bands; and (ii) this difference to different targets in a band. The latter takes according to the difference in illumination caused by the topography of a scene.

The disadvantages of the Ratio Bands are: (i) Exaggeration of noise and loss of texture due to the attenuation of shading; and (ii) Failure to distinguish between targets with similar spectral behavior. The first is caused because the shading is highly correlated in all bands, so, when shading is deleted, no loss of information concerning the topography.

Table 1 -Comparison of wavelengths (μm) of the spectral bands of sensors TM and ETM+ [START_REF] Maia | Zoneamento da Faixa Marinha do Setor Leste da Zona Costeira do Estado do Ceará[END_REF]. In FIGURE 3, the same composition, but taken by Landsat 7, these scars are also observed, but due to better image clarity, places where there is the appearance of vegetation on sites of former landslides are also distinguished. In FIGURE 4, the Ratio Bands of composition 5/7-4/3-4/1 RGB (Landsat 5), the region of Angra dos Reis, south coast of Rio de Janeiro, it is noted that the process of Ratio Bands highlights sites with exposed rocks, which are widespread in the region (A), pasture (B) and two points of recent fires (C).

CONCLUSIONS

Currently, due to countless orbiting satellites and the great technological advances, the images have been becoming better with spatial resolutions of the order of centimeters. For being older, the Landsat satellites series still present a poor spatial resolution of the order of 30 meters, making this the great disadvantage related to other satellites. However, there is a great advantage of their images being distributed for free to all individual bands. In Brazil, through the website of INPE, it is possible to make a register and download the needed images. Thereby, better costbenefit relation of the data isgenerated. Despite the high spatial resolution presented by the Landsat images, the images were good to several studies, such as landslides, mainly on the coast of Rio de Janeiro. Images of the TM sensor (Landsat 5) have a poorer clarity in relation to images obtained from ETM+ sensor (Landsat 7), which coupled with the spatial resolution of 30 meters, makes images of TM sensor less suitable for some studies. It is observed that the ETM+ sensor shows improvements in geometrical, radiometrical and clarity with respect to TM sensor. Images by Landsat 7 are more interesting for imaging with direct applications in studies of up to 1:25,000 scale. However, a disadvantage from Landsat 7 to Landsat 5 happens because the first one has no longer been in operation since 2003 and the second has been running until the present.

Figure 1 .

 1 Figure 1. Location of the study area, state of Rio de Janeiro, in the context of southeastern Brazil.

Figure 2 .

 2 Figure 2. Landsat 5 image, composition 4-7-3 (IHS-RGB Transformation) showing scars with rocky outcrops. Image width is 35km with the true north at the top.

Figure 3 .

 3 Figure 3. Landsat 7 image, composition 4-7-3 (IHS-RGB Transformation) showing scars with exposed rock (A); and grass (grass cover) (B). Image width is 35km with the true north at the top.

Figure 4 .

 4 Figure 4. Landsat 5 Image, 5/7-4/3-4/1 composition (RGB) showing scars with exposed rock (A), pasture (B), and recent local fires (C). Image width is 35km with the true north at the top.

Figure 5 .

 5 Figure 5. Landsat 7 Image, 5/7-4/3-4/1 composition (RGB) showing scarring with exposed rock (A), pasture (B) and landslides (C). Image width is 35km with the true north at the top.

ACKNOWLEDGEMENTS

The financial support of PNPD/CAPES is strongly acknowledged, for a Postdoctoral fellowship to the Associate Professor (UFRRJ) Fernando Machado de Mello at University of Sao Paulo (USP).