Guillaume Sabiron
email: guillaume.sabiron@onera.fr

Laurent Burlion
email: laurent.burlion@onera.fr

Grégory Jonniaux
email: gregory.jonniaux@astrium.eads.net

Erwan Kervendal
email: erwan.kervendal@astrium.eads.net

Eric Bornschlegl

Thibaut Raharijaona
email: thibaut.raharijaona@univ-amu.fr.

Franck Ruffier
email: franck.ruffier@univ-amu.fr.

Backup State Observer Based on Optic Flow Applied to Lunar Landing

The observer presented in this paper, which was based on the use of three minimalistic bio-inspired Visual Motion Sensors (VMS) detecting Optic Flow (OF) cues, states was intended as a backup solution in the case of Inertial Measurement Unit (IMU) failure. Contrary to most previous Guidance Navigation and Control (GNC) solutions for planetary landing, which have involved a sensor suite including an IMU, an innovative strategy is presented here for estimating states without any need for inertial measurements, based solely on information about the relative velocity of the images of the surrounding environment. A Linear Parameter Varying (LPV) observer designed on a LPV system linearized around a reference trajectory, estimates: the ventral OF, the expansion OF and the local pitch angle. A previously developed observer was applied here to a larger class of nonlinear systems by making an ingenious change of variable. Simulations performed on a lunar landing scenario yielded satisfactory performance and showed the robustness of the OF based observer to initial uncertainties and measurement noise.

I. INTRODUCTION

In most previous systems designed for the autonomous navigation of robotic systems, pose and attitude parameters have usually been measured or estimated during planetary landing [START_REF] Bryson | Vehicle model aided inertial navigation for a uav using low-cost sensors[END_REF], [START_REF] Wendel | An integrated gps/mems-imu navigation system for an autonomous helicopter[END_REF], [START_REF] Hua | Attitude estimation for accelerated vehicles using gps/ins measurements[END_REF]. However, vision based sensors and algorithms which meet the stringent weight, size and power consumption requirements of spatial applications, have recently provided new means of controlling these complex systems. Two different vision based approaches have been widely studied by performing numerical simulations:

• The first approach was based on the use of images along with information provided by sensors of other kinds or embedded knowledge of the terrain to reconstruct classical states such as velocities, attitude angles and angular velocities [START_REF] Jean-Marius | Integrated vision and navigation for a planetary lander[END_REF], [START_REF] Parkes | LIDARbased GNC for planetary landing: Simulation with PANGU[END_REF], [START_REF] Cheviron | Robust nonlinear fusion of inertial and visual data for position, velocity and attitude estimation of uav[END_REF], [START_REF] Flandin | Vision based navigation for planetary exploration[END_REF], [START_REF] Mourikis | Vision-aided inertial navigation for spacecraft entry, descent, and landing[END_REF], [START_REF] Grabe | On-board velocity estimation and closed-loop control of a quadrotor uav based on optical flow[END_REF]. Once these states have been estimated, classical control theory can be used to bring the system autonomously to the Noordwijk, The Netherlands Eric.Bornschlegl@esa.int appropriate destination. However, these recent developments are often associated with a high computational cost, mainly due to the image processing algorithm extracting visual cues from the cameras output. • In the second approach, the system was driven on the basis of relative information extracted from images of the environment. It has been established that flying insects use the Optic Flow (OF), which provides them with relative angular velocity and proximity information with respect to obstacles, to navigate swiftly in unknown complex environments. The authors of several robotic studies inspired by insects' behavior have used the OF to perform hazardous tasks such as taking off, terrainfollowing, and landing safely and efficiently [START_REF] Ruffier | Optic flow regulation: the key to aircraft automatic guidance[END_REF], [START_REF] Hérissé | A terrainfollowing control approach for a vtol unmanned aerial vehicle using average optical flow[END_REF], avoiding frontal obstacles [START_REF] Barrows | Mixed-mode VLSI optic flow sensors for in-flight control of a Micro Air Vehicle[END_REF], [START_REF] Griffiths | Maximizing miniature aerial vehicles[END_REF], [START_REF] Ruffier | Aerial robot piloted in steep relief by optic flow sensors[END_REF], [START_REF] Beyeler | OptiPilot: control of takeoff and landing using optic flow[END_REF], tracking a moving target [START_REF] Kendoul | An adaptive visionbased autopilot for mini flying machines guidance, navigation and control[END_REF] and hovering and landing on a moving platform [START_REF] Hérissé | Landing a VTOL unmanned aerial vehicle on a moving platform using optical flow[END_REF]. OF based lunar landing has been addressed in several studies using either a nonlinear observer connected to a LQ controller to track a constant OF reference in [START_REF] Valette | Biomimetic optic flow sensing applied to a lunar landing scenario[END_REF] or PID type controllers to track constant OF references [START_REF] Medici | Neuromorphic computation of optic flow data bio-inspired landing using biomorphic vision sensors[END_REF] or exponentially decreasing [START_REF] Izzo | Landing with time-to-contact and ventral optic flow estimates[END_REF] or more recently, Model Predictive Control [START_REF] Izzo | Nonlinear model predictive control applied to vision-based spacecraft landing[END_REF]. In all these studies, sensors oriented at constant angles of 90 • and sometimes 135 • were used to compute specific OF expressions [START_REF] Medici | Neuromorphic computation of optic flow data bio-inspired landing using biomorphic vision sensors[END_REF] so that:

ω x = Vx h = ω 90 • ω z = Vz h = ω 90 • -2ω 135 • (1)
Keeping (1) constant or near a slowly varying reference trajectory while h decreases ensures a soft touchdown of the closed loop system. In other studies, OF measurements have been used as a means of estimating the usual states of the system along with other more classical sensors such as Inertial Measurement Units (IMUs), sonars, Global Positioning System (GPS), airspeed sensors and/or accelerometers [START_REF] Garratt | Vision-based terrain following for an unmanned rotorcraft[END_REF], [START_REF] Kendoul | Optic flow-based vision system for autonomous 3d localization and control of small aerial vehicles[END_REF], [START_REF] Valette | Biomimetic optic flow sensing applied to a lunar landing scenario[END_REF], [START_REF] García Carrillo | Combining stereo vision and inertial navigation system for a quadrotor uav[END_REF].

However, in all of these studies, the Inertial Measurement Unit (IMU) was the crucial cue: nowadays, an IMU is the corner stone of all the autopilots designed for vehicles of all kinds. For instance, in August 2012, the NASA scientists reported the occurrence of an IMU failure which caused the very advanced Morpheus lander prototype to crash while performing its first untethered flight [START_REF] Devolites | Morpheus 1.5 a lander failure investigation results[END_REF]. IMU-less backup solutions are still urgently required in order to prevent accidents of this kind.

In the present paper, which focuses on the navigation part of the whole GNC strategy (observation issue), it is not proposed to address the important closed-loop control and guidance issues arising in Guidance Navigation and Control (GNC) design.

In this study (for the first time to the best of our knowledge), a novel backup solution was simulated for estimating the main parameters required to perform a bio-inspired planetary landing, namely the ventral and expansion OFs and the pitch angle [START_REF] Ruffier | Optic flow regulation: the key to aircraft automatic guidance[END_REF], [START_REF] Medici | Neuromorphic computation of optic flow data bio-inspired landing using biomorphic vision sensors[END_REF], [START_REF] Izzo | Landing with time-to-contact and ventral optic flow estimates[END_REF] (The expansion OF has also been referred to as the inverse of the Time To Contact) in the case of IMU failure using only three strapped down OF sensors. With these hard-mounted sensors, there is no need for gimbal systems, which usually require a dedicated IMU. No additional exteroceptive or proprioceptive sensors were used in this setup. The reference landing trajectory was used to design a Linear Parameter Varying (LPV) observer. These bio-inspired lightweight, small-sized, and energy efficient sensors featuring only 6 pixels, which were previously developed and tested outdoors onboard a UAV flying freely over various fields, gave promising results in terms of the measured OF on a real life system subjected to strong disturbances [START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF].

In section 2, the dynamic model for the lander is described and a mathematical definition of the OF is presented. Section 3 describes the fusion scheme based on the least squares of OF measurements along with a linearized version of the model around a reference trajectory. Section 4 describes how the OF based observer was designed. Section 5 presents the results of the PANGU-based simulations. Section 6 contains some concluding comments and outlines our plans for future developments.

II. LUNAR LANDER DYNAMICS AND OPTIC FLOW EQUATIONS

In this section, the dynamic model of the OF system shown in Fig. 1, and the mathematical background are described. The dynamic system studied here consisted of a spacecraft, which was actuated via the main thrust u th acting along e 3 the third canonical basis vector in the bodys fixed reference frame associated with the vector basis (e 1 ; e 2 ; e 3 ) and u θ creating a pitch torque. Since the lunar atmosphere was very thin, no friction or wind forces were applied to the lander. In line with previous studies in the literature, the lunar ground was assumed to be flat (with an infinite radius of curvature) (see [START_REF] Jean-Marius | Integrated vision and navigation for planetary exploration -final report[END_REF]). The landers dynamic motion can be described on the vertical plane (e x ; e z ) of the inertial frame I associated with the vector basis (e x ; e y ; e z ) by the following dynamic system:

           Vx (t) = sin(θ(t)) m ldr (t) u th (t) Vz (t) = cos(θ(t)) m ldr (t) u th (t) -g M oon q(t) = R I u θ (t) ṁldr (t) = -u th (t) Isp th .g Earth + -|u θ (t)| Isp θ .g Earth (2)
where V x,z denotes the lander's velocities in I, m ldr stand for the lander's mass, which was assumed to be measured at all times, θ is the pitch angle (the pitch angular rate is written q = dθ dt since this is taken to be a simplified 2-D problem), t θ e 1 e 3 u th

FoE Φ N Φ 1 ω α N +θ ω α 1 +θ ω α 2 +θ D N D 1 D 2 γ α 1,2,••• ,N V x V V z Inertial Frame e x
e z e y Fig. 1. Definition of the body-fixed frame B, the inertial frame I and notations for the lander's dynamics and the optic flow. Example of a N VMS setup. Sensors fixed to the landers structure measured the OF ω α N +θ depending on the lander's position, velocity, pitch angle and pitch angular velocity (see [START_REF] Bryson | Vehicle model aided inertial navigation for a uav using low-cost sensors[END_REF]). It can be noted that Φ(t) = α+θ(t)+γ(t). The canonical vector e 2 of B was not included here for the sake of clarity.

denotes the time, and g M oon denotes the lunar acceleration due to the gravity (g M oon = 1.63 m/s 2 , g M oon was taken to be constant due to the low initial altitude). I is the landers moment of inertia, and R denotes the eccentricity of the attitude thrusters from the center of mass. I sp is the specific impulse: I sp th = 325s that of the braking thrusters, I sp θ = 287s that of the attitude thrusters and g Earth = 9.81 m/s 2 is the Earth's gravity. In the first step, we focused on a 2-D setup (planar motion) in line with previous studies and because of the 6-pixel OF sensors capabilities presented in [START_REF] Sabiron | Low-speed optic-flow sensor onboard an unmanned helicopter flying outside over fields[END_REF]. In order to progress to a 3-D setup, we would have to adopt a larger photoreceptor chip with a matrix-shaped design to be able to measure the 2-D OF involved in roll and yaw movements.

In the vertical plane, the OF ω(Φ), the sole value measured in this study, was defined in line with [START_REF] Koenderink | Facts on optic flow[END_REF] as follows:

ω(Φ) = V D sin(Φ) -q (3)
where the term V D sin(Φ), which has been called the translational OF, depends on the linear velocity V expressed in the inertial frame, the distance from the ground D in the gaze direction and the elevation angle Φ (i.e. the angle between the gaze direction and the heading direction). In order to use the useful properties of the translational OF, a bulky gyro stabilized gimbal system has often been used to compensate for spacecrafts rotations and thus to keep the visual sensors oriented in the same direction. This means not only that the last term in (3) corresponding to the rotational OF no longer occurs in the measurements, but also that the angle Φ -γ was kept constant which simplifies calculations.

Although this is simple, useful mathematical framework no gimbal systems were used in this paper, it would have required an extra IMU, which is not possible since we are dealing with a case of IMU failure. The challenge of using the present strapped down sensor setup was twofold.

First the angle Φ -γ was no longer constant, which ruled out simple calculations such as those performed in [START_REF] Argyros | Biomimetic centering behavior [mobile robots with panoramic sensors[END_REF]. Secondly, the rotational OF ω R = -q was included in the measurements. This component ω R was usually subtracted from the measured OF, ω meas , using IMU measurements: this operation is known as the derotation process (see [START_REF] Argyros | Biomimetic centering behavior [mobile robots with panoramic sensors[END_REF]). Unfortunately, since the IMU was lost, it was no longer possible to perform the derotation process or to close the attitude control loop depending on pitch and angular pitch velocity measurements. We therefore used an OF based observer to estimate ω x and ω z during the descent as well as the pitch angle θ without using any gimbaled sensors, IMUs or velocity measurements. In the case of IMU failure, this solution could be used as a backup solution to feed control laws with accurate estimates and enable a small airborne vehicle to perform a soft landing.

III. OPTIC FLOW FUSION AND LPV MODEL DEFINITION

A fusion scheme was applied to the OF measurements and the landers dynamics (2) were linearized around a reference trajectory so as to obtain a Linear Parameter Varying (LPV) state space model, which was used throughout to design the observer.

Let us take only N available measurements (ω α1 ; • • • ; ω α N ) to estimate the pitch angle θ and useful values of OF, namely the ventral OF, ω x , and the expansion OF, ω z :

ω x (t) ω z (t) = Vx h(t) Vz h(t)
from the available measurements.

Assuming that we are dealing with a practically flat ground (i.e. D = h/ cos( π 2 -Φ + γ), where γ denotes the flight path angle (the orientation of the velocity vector with respect to the local horizontal), h denotes the ground height, and Φ -γ denotes the angle between the gaze direction and the local horizontal) and using the notations presented on Fig. 1, we obtain the following general expression for the 2D OF:

ω Φ (t) = V (t) D(t) sin(Φ(t)) -q(t) (4)
where

Φ(t) = α + θ(t) + γ(t), and D = h(t) sin(α+θ(t))
Based on (1) this gives:

ω α+θ (t) = 1 2   1 -c (2 (α + θ(t))) s (2 (α + θ(t))) -2   T   ω x (t) ω z (t) q(t)   (5)
where sin(γ(t)) = Vz V and cos(γ(t)) = Vx V and α, the only time invariant parameter (time notations were dropped for the sake of clarity); s (•) and c (•) denotes the sine function and the cosine function.

Extracting a time invariant matrix depending on α yields

ω α+θ =   1 2 c (2α) s (2α)   T   ωx -2q -c(2θ)ωx + s (2θ) ωz s(2θ)ωx + c (2θ) ωz   (6)
The following linear system can then be solved using N ≥ 3 measurements if H T H is invertible:

  ω x -2q -cos(2θ)ω x + sin (2θ) ω z sin(2θ)ω x + cos (2θ) ω z   = H -1 lef t N    ω α1+θ . . . ω α N +θ    (7) where the left inverse H -1 lef t is defined as H -1 lef t = (H T H) -1 H T , with H -1 lef t H = I n , H ∈ M (m, n, R
), m > n and H defined as follows:

H N =    1 cos (2α 1 ) sin (2α 1 ) . . . . . . . . . 1 cos (2α N ) sin (2α N )   
Remark It should be noted that for N=3 the matrix H matrix is a square matrix:

  ω x -2q -c(2θ)ω x + s (2θ) ω z s(2θ)ω x + c (2θ) ω z   = H -1   ω α1+θ ω α2+θ ω α3+θ   ( 8
)
with

H =   1 2
cos (2α 1 ) sin (2α 1 )

1 2 cos (2α 2 ) sin (2α 2 ) 1 2 cos (2α 3 ) sin (2α 3 )  
It has to be checked that the following condition is satisfied to ensure that H is invertible

s (2 (α 3 -α 2 )) + s (2 (α 1 -α 3 )) + s (2 (α 2 -α 1 )) = 0 (9)
The result of the linear least squares calculation with N > 3 (and the matrix inversion with N = 3) gives a nonlinear system with three equations and four unknowns (i.e. ω x , ω z , θ and q) which it is impossible to solve analytically. We therefore had to make use of the knowledge available about the systems dynamics, which was done by designing an LPV observer for a linearized model of the lander's dynamics. Roughly speaking, we therefore linked together two unknowns θ and q, since θ = q. During space missions including entry, descent and landing phases, a reference trajectory has often been provided. This trajectory has to be followed during the actual landing to ensure the safety of the lander and reduce the fuel consumption (see [START_REF] Liu | Robust guidance and control of lunar lander using model reference approach[END_REF], [START_REF] Zhou | Optimal guidance for lunar module soft landing[END_REF]). The reference trajectory can be expressed in terms of state trajectories, crater patterns, or a 2-D or 3-D image database [START_REF] Johnson | Overview of terrain relative navigation approaches for precise lunar landing[END_REF]. In this study, we adopted a landing scenario with a computed offline reference trajectory which had to be followed. Our main assumption was that we would stay sufficiently close to this reference trajectory to be able to derive and use an LPV model for the system around the reference trajectory defined by X * = [h * , ω *

x , ω * z , θ * , q * ] T involving the input control sequences u * = [u * th , u * θ ] T from the system dynamics (2) and the results of the linear least squares calculations [START_REF] Flandin | Vision based navigation for planetary exploration[END_REF]. Linearized outputs are given by

Y =   1 0 0 -c (2θ * ) s (2θ * ) 2ω * x s (2θ * ) + 2ω * z c (2θ * ) s (2θ * ) c (2θ * ) 2ω * x c (2θ * ) -2ω * z s (2θ * )     ∆ωx ∆ωz ∆θ   +   -2 0 0   ∆q +   ω * x -2q * -c (2θ * ) ω * x + s (2θ * ) ω * z s (2θ * ) ω * x + c (2θ * ) ω * z   (10)
Which was rewritten as follows

Y = C 1 ∆h ∆q + C 2 (X * (t))   ∆ω x ∆ω z ∆θ   + Y * (11)
where

C 1 =   0 -2 0 0 0 0   .
A linearized state space model around a reference trajectory is given by

             ∆ Ẋ =      ∆ ḣ ∆ q ∆ ωx ∆ ωz ∆ θ      = A (ρ (t)) ∆X + B (ρ (t)) ∆u th ∆u θ ∆Y = C (ρ (t)) ∆X = Y -Y * (12)
with

A =         ω * z 0 0 h * 0 0 0 0 0 0 -s(θ * )u * th m ldr h * 2 0 -ω * z -ω * x c(θ * )u * th m ldr h * -c(θ * )u * th m ldr h * 2 + g M oon h * 2 0 0 -2ω * z -s(θ * )u * th m ldr h * 0 1 0 0 0         (13)
B =        0 0 0 R/I s(θ * ) m ldr h * 0 c(θ * ) m ldr h * 0 0 0        ; C (ρ (t)) = C1 C2 (ρ (t)) (14)
The time varying vector ρ(t) depended on the reference trajectory, its associated input control signals and on the lander's mass so that ρ(t) = [h * , ω *

x , ω * z , θ * , q * , u * th , m ldr ] T . An LPV system was obtained; an LPV observer will now be designed to estimate state deviations from the reference trajectory, assuming that these deviations are small.

IV. OBSERVER DESIGN FOR A CLASS OF LPV SYSTEMS AND APPLICATION TO A LUNAR LANDING SCENARIO

In this section, we present an LPV observer based on the solution proposed in [START_REF] Besanc ¸on | Observer synthesis for a class of nonlinear control systems[END_REF], which was extended to include a larger class of nonlinear systems and applied to the lunar landing scenario adopted in this paper.

A. LPV Observer Design

Let us now consider an LPV system having the following form

Ẋ (t) = A (ρ (t)) X + B (ρ (t)) u (t) y (t) = C (ρ (t)) X (t) = C 1 C 2 (ρ (t)) X (t) (15)
where

X ∈ R n , ρ ∈ R m , u ∈ R p , y ∈ R r , A : R m → M(n, n), B : R m → M(n, p), C : R m → M(r, n), C 1 ∈ M(r, n -s), C 2 :
R m → M(r, s) where M(k, l) denotes the space consisting of k × l matrices with k rows and l columns with coefficients provided in R. I a denotes the identity matrix having the size a × a.

It is worth noting that the system described by ( 11)-( 12) belongs to this class of LPV systems.

Lemma 1 Let us consider the system [START_REF] Hérissé | Landing a VTOL unmanned aerial vehicle on a moving platform using optical flow[END_REF] and assume that the parameter ρ(t) is measured and regularly persistent (see [START_REF] Besanc ¸on | Observer synthesis for a class of nonlinear control systems[END_REF] for a definition) and the sub-matrix C 2 (ρ(t)) is invertible and differentiable with respect to time, then for all gain matrices Θ > 0,the system:

Ẋ (t) = Ā X (t) + Bu (t) -S -1 CT (ŷ (t) -y (t)) (16) Ṡ (t) = -ΘS (t) -Ā (ρ(t)) T S (t) -S (t) Ā (ρ(t)) + C (ρ(t)) T C (ρ(t)) S(0) > 0 (17)
where

X (t) = I n-s 0 n-s,s 0 s,n-s C 2 (ρ(t)) X (t) = M (ρ(t)) X (t)
which yields

Ā = Ṁ (ρ(t)) M (ρ(t)) -1 + M (ρ(t)) A (ρ(t)) M (ρ(t)) -1 B = M (ρ(t)) B (ρ(t)) ; C = C 1 I s
is an observer for [START_REF] Griffiths | Maximizing miniature aerial vehicles[END_REF].

Proof: Let us now consider the new state vector resulting from the following change of variable X = M (ρ(t)) X: [START_REF] Griffiths | Maximizing miniature aerial vehicles[END_REF] becomes

Ẋ = Ṁ (ρ(t)) M (ρ(t)) -1 + M (ρ(t)) A (ρ(t)) M (ρ(t)) -1 X + M (ρ(t)) B (ρ(t)) u and y = C X
which corresponds to the class of systems covered by theorem 2.1 presented in [START_REF] Besanc ¸on | Observer synthesis for a class of nonlinear control systems[END_REF] (additional proof can be found in [START_REF] Hammouri | Observer synthesis for stateaffine systems[END_REF]). As previous authors have pointed out, with this observer, we ensure that ǫ (t)

2 ≤ ae -λmin(Θ)t where ǫ (t) = X (t)-X (t) is the estimation error, λ min (Θ) is the smallest eigenvalue of Θ and a is a constant depending on the initial errors, Θ and u.

B. Application to the lunar lander LPV state space model

Applying the change of variable to the system described in ( 11)-( 12) gives the new state vector

∆ X =       ∆h ∆q ∆x 1 ∆x 2 ∆x 3       = I 2 0 2,3 0 3,2 C 2 (ρ (t)) ∆X ∆ X = M (ρ (t)) ∆X
where 0 a,b denotes the null matrix having the size a × b.

The measurement equation can then be written with a time invariant observation matrix C

Y = Y * + C 1 I 3 ∆ X = Y * + C∆ X
and the state equation can be written as follows

∆ Ẋ = I 2 0 2,3 0 3,2 Ċ2 C -1 2 + C 2 AC -1 2 ∆ X + I 2 0 2,3 0 3,2 C 2 B∆u with Ċ2 =   0 0 0 2q * s (2θ * ) 2q * c (2θ * ) C 223 2q * c (2θ * ) -2q * s (2θ * ) C 233   C2 23 = 2 2q * c 2θ * ω * x + s 2θ * ω * x -2q * s 2θ * ω * z + c 2θ * ω * z C2 33 = 2 -2q * s 2θ * ω * x + c 2θ * ω * x -2q * c 2θ * ω * z -s 2θ * ω * z
By applying lemma (1), the observer can be expressed as follows:

∆ Ẋ = I2 02,3 03,2 Ċ2C -1 2 + C2AC -1 2 ∆ X + I2 02,3 03,2 C2 B∆u -S -1 C T Ŷ -Y = Ā X (t) + Bu (t) -S -1 CT Ŷ (t) -Y (t) (18)
Ṡ = -ΘS -ĀT S -S Ā + CT C; S(0) > 0 (19)
Θ does influence the observers convergence time, but a high gain Θ is liable to amplify any measurement noise present in the simulation involving virtual images of the lunar ground. We therefore used a trial and error method and engineering knowledge of the dynamics and sensors to define the S and Θ matrices as follows:

S(0) =    1 0.2 I 2 0.2    and Θ =    20 0.05 20I 2 0.3   
Θ is a tuning parameter which affects the rate of convergence of the estimated states.

V. SIMULATION RESULTS WITH 3 OF SENSORS (N = 3)

The results of the simulation illustrate the performances of the LPV observer. We focused here on performing simulations based on real image processing rather than using methods based on the intentional corruption of signals by adding Gaussian noise. Simulations were run under open loop conditions to obtain realistic OF measurements using three sensors placed on the lander at random angles α 1 = 90 • , α 2 = 120 • and α 3 = 160 • in a range corresponding to the landing scenario adopted here, so that they remained ground oriented during the entire descent phase and the condition defined by ( 9) was satisfied. These small inter-receptor and acceptance angles make it possible to compute very low velocities. Since we have such a narrow field of view, even high spatial frequency contrasts will be detected by the photodiodes, which is extremely useful at low OF levels, where fewer contrasts occur in the sensor's line of sight. In the simulated VMS model, the photoreceptors' output was simulated at each time step by convolving the PANGU-generated lunar surface image with the 2-D Gaussian filter. The simulated 6-pixel VMSs based on the actual code implemented in the sensor then assessed the OF. Fig. 2 shows a simulated image of the lunar ground generated using PANGU. One can see various craters, boulders and shadows caused by the rims of craters and the elevation of the sun. Simulations were run in open loop with precomputed input control signals using a lunar landing scenario starting with 3. Evolution of estimated ωx, ωx and θ in the nonlinear system with the newly designed LPV observer. The three sub-figures show the actual state of the nonlinear system (blue line), the reference state (dashed black line), the open loop linearized system (△-mark line) and the estimated states ( * -mark line). The open loop linearized system states correspond to a simulation where the linearized system is run in open loop from the same initial conditions as the observer (states are denoted ωx OLlin , ωz OLlin and θ OLlin ). One can see that the estimated states converged quickly toward the actual states of the nonlinear system, which were deliberately intended not to follow the reference trajectory. It should be noted in addition that the open loop linearized system alone would not have sufficed to obtain an accurate state estimation. Estimated states were initialized, giving ∆ωx 0 = -0.126 • /s, ∆ωz 0 = 0.23 • /s, ∆ θ0 = 3.23 • .

h 0 = 1496m, ω x0 = 2.18 • /s, ω z0 = Time (s) c) θ(t f ) -θ OLlin (t f ) =-3.7 • /s θ(t f ) -θ(t f ) =1.1 • /s θref (--), θ, θ ( * ), θ OLlin (△), ( • ) b) ω z (t f ) -ω zOLlin (t f ) =-14.6 • /s ω z (t f ) -ωz (t f ) =0.3 • /s ω zref (--), ω z , ωz ( * ), ω zOLlin (△), ( • /s) a) ω x (t f ) -ω xOLlin (t f ) =3.3 • /s ω x (t f ) -ωx (t f ) =-0.6 • /s ω xref (--),
-1.72 • /s, θ 0 = -53.43 • , q 0 = -0.92 • /s, m ldr0 = 758kg. Initial states deviations from the reference trajectory were chosen so that ∆h 0 = 100m, ∆ω x0 = -0.115 • /s, ∆ω z0 = -0.08/s, ∆θ 0 = -2 • , ∆q 0 = -0.2 • /s. The estimated states were initialized, giving ∆ ĥ0 = 60m, ∆ω x0 = -0.126 • /s, ∆ω z0 = 0.23/s, ∆ θ0 = 3.23 • , ∆q 0 = -0.2 • /s. Fig. 3 and Fig. 4 present simulations based on the state estimation technique presented here. In sub-figures 3.a-c, one can see the actual state (blue line), the reference state (dashed black line), the open loop linearized system (△-mark line) and the estimated states ( * -mark line). In the sub-figures 4.c, one can see the actual outputs (blue line), the reference outputs (dashed black line), the linearized open loop outputs (△-mark line), the estimated outputs ( * -mark line) and the output given by the PANGU-based measurements (▽-mark line). The open loop linearized system states and outputs were obtained in simulation with the linearized system run in open loop from the same initial conditions as the observer (states are denoted ω x OLlin , ω z OLlin , θ OLlin , h OLlin , q OLlin and y OLlin ).

As expected ωx , ωz and θ converged quickly toward the actual states, and it is worth noting that we also obtained a raw estimates of h which took longer to converge but were not intended to be used in the future control scheme. q seems to have been more sensitive to the simulated imagebased measurements than the other estimates: this was due to the events created by the OF measurement techniques . Evolution of ĥ, q and the outputs of the observer Ŷ = [ŷ 1 , ŷ2 , ŷ3 ] T in the nonlinear system equipped with an LPV observer. One can see from this figure that the estimated states gave good raw estimates of the actual states despite the intrinsic sensitivity of the system to noise during the last few seconds. The estimated outputs were very close to the actual outputs of the nonlinear system. The estimated states were initialized in such a way that ∆ ĥ0 = 60m, ∆q 0 = -0.2 • /s. The open loop linearized system states correspond to a simulation where the linearized system is run in open loop from the same initial conditions as the observer (states are denoted h OLlin , q OLlin and y OLlin ).

y 3 (t f ) -ŷ3 (t f ) =0.6 • /s, y 2 (t f ) -ŷ2 (t f ) =0.3 • /s, y 1 (t f ) -ŷ1 (t f ) =0.2 • /s, Time (s) y ref , y, ŷ ( * ), y OLlin (△), y meas (▽), ( • /s) b) q(t f ) -q OLlin (t f ) =-0.4 • /s q(t f ) -q(t f ) =-0.4 • /s q ref (--), q, q ( * ), q OLlin (△), ( • /s) a) h(t f ) -h OLlin (t f ) =-470.7m h(t f ) -ĥ(t f ) =-68.3m h ref (--), h, ĥ (
as well as the noise in the measurements. Simulations with theoretical measurements were also run to check the validity of this statement (i.e. q -q tended to 0 when the outputs were perfect). Open loop linearized system (△-mark line) states (without the presence of an observer) did not converge toward true values, and even diverged with ω x and ω z .

VI. CONCLUSIONS

This study addressed the state estimation problem by using OF sensors without any need for IMU measurements. In the novel LPV observer adapted from [START_REF] Besanc ¸on | Observer synthesis for a class of nonlinear control systems[END_REF] to our class of nonlinear systems by making a change of variable, the measurements used required only three lightweight bioinspired visual motion sensors hard-mounted onto the landers structure (i.e. ω α1+θ , ω α2+θ and ω α3+θ ). No inertial measurements (attitude, angular velocities, angular or linear accelerations) or linear velocities or even altitudes were needed in this setup. The promising results obtained here showed the effectiveness of the observer in simulations based on software-generated images of the lunar ground. The ventral OF, expansion OF and pitch angle estimates were very near the actual states although substantial sensor noise and the estimated initial states were slightly inaccurate. Due to the extreme minimalism of the present OF sensors and the overall backup solution, the performances of the present observer are not comparable with the high-accuracy IMU devices available. However, the present observer was intended as a backup means of driving the lander safely toward the lunar surface based on the use of small, energy efficient sensors if major sensor failures of any kind (involving not only the IMU) should occur.

The next step will consist in designing the control laws, taking ω x , ω z and θ as inputs to bring the lander safely to the ground while following a precomputed reference trajectory.

To overcome the poor estimation of the pitch angular speed, two approaches could be used:

• controlling the lander's attitude only using pitch measurement (basically a position feedback), • or combining the estimated angular pitch rate with the time derivative of the estimated pitch angle. It would also be interesting to compare results obtained with the present LPV observer with traditional extended unscented Kalman filters, which are designed on very similar lines and are being increasingly used to meet challenges of this kind. The LPV observer was chosen because of its wellestablished success with nonlinear systems, and also because of the proof of convergence provided by the Lyapunov theory underlying Besanc ¸on, Bornard and Hammouri's observer [START_REF] Besanc ¸on | Observer synthesis for a class of nonlinear control systems[END_REF]. A cost-minimization scheme could be developed in order to improve the settings of the observers matrices Θ and S(0). Further simulations will be performed in order to define the limitations of the solution proposed in terms of the initial errors and measurement noise.

Fig. 2 .

 2 Fig. 2. Example of an image of the lunar ground obtained using PANGU software 2.70 PANGU software was used to generate images of the lunar surface, taking the position of the system, the elevation of the sun and the camera's properties into account. The simulated lunar surface was irregular and sometimes included craters up to 40 m deep. The images generated by PANGU contained 256 gray-scale levels and had a resolution of 256×256 pixels. Each of the OF sensors (also called VMSs) included six photoreceptors: the visual axes of each pair of photoreceptors were separated by the inter-receptor angle ∆ϕ = 0.1 • . The angular sensitivity of each photoreceptor obeyed a 2-D Gaussian function mimicking the angular sensitivity of the fly's photoreceptors with the acceptance angle (the angular width at half height) ∆ρ = ∆ϕ = 0.1 • .These small inter-receptor and acceptance angles make it possible to compute very low velocities. Since we have such a narrow field of view, even high spatial frequency contrasts will be detected by the photodiodes, which is extremely useful at low OF levels, where fewer contrasts occur in the sensor's line of sight. In the simulated VMS model, the photoreceptors' output was simulated at each time step by convolving the PANGU-generated lunar surface image with the 2-D Gaussian filter. The simulated 6-pixel VMSs based on the actual code implemented in the sensor then assessed the OF. Fig.2shows a simulated image of the lunar ground generated using PANGU. One can see various craters, boulders and shadows caused by the rims of craters and the elevation of the sun. Simulations were run in open loop with precomputed input control signals using a lunar landing scenario starting with h 0 = 1496m, ω x0 = 2.18 • /s, ω z0 =

  Fig.3. Evolution of estimated ωx, ωx and θ in the nonlinear system with the newly designed LPV observer. The three sub-figures show the actual state of the nonlinear system (blue line), the reference state (dashed black line), the open loop linearized system (△-mark line) and the estimated states ( * -mark line). The open loop linearized system states correspond to a simulation where the linearized system is run in open loop from the same initial conditions as the observer (states are denoted ωx OLlin , ωz OLlin and θ OLlin ). One can see that the estimated states converged quickly toward the actual states of the nonlinear system, which were deliberately intended not to follow the reference trajectory. It should be noted in addition that the open loop linearized system alone would not have sufficed to obtain an accurate state estimation. Estimated states were initialized, giving ∆ωx 0 = -0.126 • /s, ∆ωz 0 = 0.23 • /s, ∆ θ0 = 3.23 • .

Fig. 4

 4 Fig.[START_REF] Beyeler | OptiPilot: control of takeoff and landing using optic flow[END_REF]. Evolution of ĥ, q and the outputs of the observer Ŷ = [ŷ 1 , ŷ2 , ŷ3 ] T in the nonlinear system equipped with an LPV observer. One can see from this figure that the estimated states gave good raw estimates of the actual states despite the intrinsic sensitivity of the system to noise during the last few seconds. The estimated outputs were very close to the actual outputs of the nonlinear system. The estimated states were initialized in such a way that ∆ ĥ0 = 60m, ∆q 0 = -0.2 • /s. The open loop linearized system states correspond to a simulation where the linearized system is run in open loop from the same initial conditions as the observer (states are denoted h OLlin , q OLlin and y OLlin ).

ACKNOWLEDGEMENTS

We thank J. Blanc for improving the English manuscript. The authors would like to thank the three anonymous referees for their valuable comments and suggestions to improve the quality of the paper.

* This research work is co-funded by CNRS Institutes (Life Science; Information Science; Engineering Science and Technology), the Aix-Marseille University, European Space Agency, ONERA the French Aerospace Lab and Astrium Satellites under ESA's Networking/Partnering Initiative program (NPI) for advanced technologies for space.


S Calloch

K Taillard

S Arbab Chirani

C Lexcellent

E Patoor

Relation between the martensite volume fraction and the equivalent transformation strain in shape memory alloys

Keywords: Superelasticity, Multiaxial Testing, Resistivity, model-experiment dialogue

This study concerns the pseudoelasticity of SMA. A series of tests under multiaxial loadings and a micro-macro model are used to show the validity of a conjecture concerning the relation between the volume fraction of martensite and the equivalent transformation strain.

INTRODUCTION

Superelasticity is a particular behavior of Shape Memory Alloys (SMA) [START_REF] Patoor | Les alliages à mémoire de forme[END_REF]. It is observed when SMA is mechanically tested at a temperature greater than the A f (austenite finish) temperature. In this case, SMA show reversibly hysteretic behavior up to strain levels near 8%. Recently, a macroscopic model which permits to simulate the superelasticity of SMA under complex multiaxial loadings has been proposed [START_REF] Bouvet | A phenomenological model for pseudoelasticity of SMA under multiaxial proportional and nonproportional loadings[END_REF]. The originality of this model is that, the elastic domain of the material during its two phased state is given by the intersection of two transformation surfaces. The first one permits to control the direct transformation (i.e., austenite _ martensite), and the second one the inverse transformation (i.e., martensite _ austenite). The reorientation mechanism, which is produced just under nonproportional multiaxial loading, is activated when the loading point is situated at the intersection of these two transformation surfaces. None supplementary variable is introduced in the model, the reorientation is modelized (at the adopted description level) by a simultaneous direct and inverse transformation. This model has been validated by a large base of multiaxial testing. Always in this model, a conjecture has been proposed concerning the relation between the equivalent transformation strain and the martensite volume fraction. Up to now, this conjecture has been validated just in the case of an uniaxial tension loading [START_REF] Vacher | Study of pseudoelastic behavior of polycrystalline SMA by resistivity measurements and acoustic emission[END_REF]. The objective of the study presented in this paper is to realize a series of tension-torsion proportional testing with the aim of discussing about the validity of the conjecture. The first part of the paper concerns the description of the realized experiments and the obtained results. We can see that a linear relation between the martensite volume fraction and the equivalent transformation strain is really obtained. The proportionality coefficient seems depending on the loading direction. With the aim of analyzing these results, a micromechanical model developed by Patoor et al. [START_REF] Patoor | Micromechanical modeling of thermoelastic behavior of shape memory alloys[END_REF][START_REF] Entemeyer | Etude micromécanique du comportement thermomécanique des alliages à mémoire de forme[END_REF], is used. This more physics model operating as a virtual testing machine permits to discuss the validity of the conjecture.

EXPERIMENTAL RESULTS UNDER MULTIAXIAL LOADING

The studied material is a polycrystal of Cu-Al-Be SMA. Its composition is Cu 87.75 wt%, Al 11.33 wt%, Be 0.45 wt%. The tension-torsion tubular samples have been machined in the drawn bars then heat treated as the following: heated during 20 minutes at 650 °C then quenched in bowling water during 1 hour. The four characteristic transformation temperatures have been determined using the resistivity measurement. They are A s = 19°C, A f = 32°C, M s = 21°C and M f = -7°C.

EXPERIMENTAL DETERMINATION OF THE TRANSFORMATION YIELD SURFACE

A series of experiments under tension (compression)-internal pressure, tension (compression)-torsion and bi-compression loading have been realized with the aim of determining the form of the transformation yield surface of the studied alloy. Detailed information about this study can be found in Bouvet et al. [START_REF] Bouvet | Mechanical behavior of a Cu-Al-Be SMA under multiaxial proportional and nonproportional loadings[END_REF]. Just the necessary results are reported in this paper. The figures 1 and 2 show one part of these experimental results.The results of these experiments permitted to define equivalent stresses for isotropic and anisotropic SMA [START_REF] Bouvet | Mechanical behavior of a Cu-Al-Be SMA under multiaxial proportional and nonproportional loadings[END_REF]. In the case of an isotropic material, the following equivalent stress has been defined: where σ and y σ are respectively the equivalent von Mises stress and the third invariant of the stress tensor defined as below:

σ eq iso = σ eq iso σ ,y σ ( ) = σ .g y σ ( ) (1)
σ = 3 2 σ D : σ D et y σ = 27 2 det σ D ( ) σ 3 (2)
and g y σ ( ) is a function defined by relation (3):

g y σ ( ) = cos 1 3 .cos -1 1-a. 1-y σ ( ) ( )       ( 3
)
where a is a material parameter which can be identified by a compression test.

An equivalent transformation strain can be defined from the equivalent stress definition, the notion of equivalence in the transformation power ( P eq tr = σ eq .˙ ε eq tr = σ : ˙ ε tr = P tr ) and the normality hypothesis

( ˙ ε tr = ˙ λ ∂σ eq ∂σ ):
ε eq tr = ε tr . g -y ε tr ( )

g -1 ( ) (4)
where ε tr is von Mises equivalent strain ε tr = 2 3 ε tr :

ε tr    
  and y ε tr is the third invariant of the transformation strain tensor defined as following:

y ε tr = 4 det ε tr ( ) ε tr ( ) 3 (5)

RELATION BETWEEN THE VOLUME FRACTION OF MARTENSITE AND THE EQUIVALENT TRANSFORMATION STRAIN

A great number of authors have proposed a linear relation between the martensite volume fraction, z, and the transformation strain in the case of the uniaxial tension [START_REF] Vacher | Study of pseudoelastic behavior of polycrystalline SMA by resistivity measurements and acoustic emission[END_REF], [START_REF] Sun | Micromechanics modelling for the constitutive behavior of polycrystalline shape memory alloys[END_REF]. We propose to generalize this relation to the all multiaxial proportional loading paths with the help of the relation ( 6):

ε eq tr = γ.z (6) where γ is a material parameter. A series of tests under tension (compression)-torsion loading using tubular samples have been realized to verify the relation [START_REF] Entemeyer | Etude micromécanique du comportement thermomécanique des alliages à mémoire de forme[END_REF]. During each proportional loading path, the martensite volume fraction and the equivalent transformation strain have been determined respectively by the resistivity and the strain measurement. The figure 3 shows the instrumentation of the tension-torsion samples schematically. The figure 4 shows the volume fraction evolution with the equivalent transformation strain measured during the 8 proportional tests realized in ε 11 , γ 3

(

) plane. On this figure it can be mentioned that until the tested level of strain, the martensite volume fraction seems to vary linearly with the equivalent transformation strain. On the other hand, the slope of this linear relation seems depending on the loading path. With the aim of analyzing these results, the Patoor et al. [START_REF] Patoor | Micromechanical modeling of thermoelastic behavior of shape memory alloys[END_REF] micromechanical model will be used in the continuation of this work.

MICROMECHANICAL MODEL

It is based on the self consistent scale transition method [START_REF] Entemeyer | Etude micromécanique du comportement thermomécanique des alliages à mémoire de forme[END_REF]. The representative elementary volume is defined at the polycrysal level. Loading is calculated at grain scale using the localization tensor.

Considering the local loading state, the martensite activable variants are detected and the volume fraction of each variant and also the induced transformation strain are evaluated. This operation is realized for all the grains of the polycrystal. After determination of the behavior of all the grains, the macroscopic behavior of the SMA is calculated using the average operations.

The data of this model are the crystallographic texture, the transformation temperatures, the material coefficients identified from the state diagram of the SMA, the normal of habit planes and the transformation directions of the 24 variants. At least an interaction matrix which defines two types of interaction (weak and strong) between the formed martensite variants in a grain. All these parameters are physical. They do not depend on the past of the material, except the crystallographic texture.

The crystallographic texture describing the different grains orientation permits to take into account the isotropy or anisotropy of the material related to the forming process [START_REF] Chirani | Superelastic behaviour modelling in shape memory alloys[END_REF]. It is well known that the initial crystallographic texture is an important parameter in the behavior of SMA: Bhattacharya [START_REF] Bhattacharya | Symmetry, Texture and the Recoverable Strain of Shape Memory Polycrystals[END_REF], Gall [START_REF] Gall | The role of texture in tension-compression asymmetry in polycrystalline NiTi[END_REF].

Using the micromechanical model permits to follow the microstructure evolution during loading and to have an idea about the martensite volume fraction level at different scale.

The studied material is a CuAlBe alloy with the macroscopic characteristics described before. In these alloys which are more stable comparing to the other copper based alloys, there are 24 martensite variants forming 6 auto accommodating groups. They are characterized by the normal to the habit plane and the transformation direction. They are the [1 1 4] type [START_REF] Moreau | LETAM[END_REF]. The transformation strain amplitude, g, is the same for all the variants and it is about 0.23. The parameters of the interaction matrix are: 50 MPa describing the weak interaction and 260 MPa for the strong one [START_REF] Patoor | Micromechanical modeling of thermoelastic behavior of shape memory alloys[END_REF]. The appearance of the martensitic transformation is related to the chemical potential ΔG which has a linear form around the M s temperature:

ΔG ch = B(T-M s ) (7)
The parameter B value is 0.23 MPa/K in equation [START_REF] Berveiller | Modelling of the plastic behavior of inhomogeneous media[END_REF] The elasticity is supposed to be isotropic and the same in the two phases.

In the following part, the micromechanical model is used with the aim of validating the shape of the transformation surfaces obtained by the macroscopic model and discussing the validity of the relation between the martensite volume fraction and the equivalent transformation strain. All the simulations, realized using the micromechanical model, are obtained from an isotropic polycrystal containing 1000 grains.

TRANSFORMATION SURFACES

The figure 5 and6 show the comparison between the transformation surfaces in the σ 11 ,σ 22 ( ) and σ 11 , τ ( ) planes calculated by the micromechanical and the macroscopic model. A good agreement between the two models is observed.

RELATION BETWEEN THE MARTENSITE VOLUME FRACTION AND THE EQUIVALENT TRANSFORMATION STRAIN

The other interesting aspect of the realized simulations concerns the equivalent stress evolution as a function of the equivalent transformation strain and the martensite volume fraction evolution as a function the equivalent transformation strain. The figure 7 and8 show these evolutions for 8 different proportional loading paths. We can consider that for the all selected loading, the curves giving the equivalent stress as a function of transformation strain are the same. This permits the validation of the defined equivalent stress and transformation strain in the case of isotropic SMA. We can also observe, except for the transformation beginning, that the martensite volume fraction evolution with the equivalent transformation strain is linear. In the case of an isotropic SMA the proportional coefficient does not depend on the considered loading path. The existence of a different slope at the beginning of the transformation had already been observed by Entemeyer [START_REF] Entemeyer | Etude micromécanique du comportement thermomécanique des alliages à mémoire de forme[END_REF]. It persists, in the material behavior, until that the transformation phase concerns the majority of the grains.

4.

CONCLUSION A testing campaign under tension-torsion loading has been developed with the aim of validation the relation between the martensite volume fraction and the equivalent transformation strain so frequently used in the macroscopic models. The results show that this relation is perfectly linear but the slope seems depending on the loading path. Using micromechanical model of Patoor et al. [START_REF] Patoor | Micromechanical modeling of thermoelastic behavior of shape memory alloys[END_REF] has permitted the study of these experimental results differently. By this way the equivalent stress and transformation strain definitions have been validated. The conjecture in the case of isotropic SMA has also been verified. The experimental results interpretation is still remained. Taking into account the anistropy of the studied alloy will maybe allow the generalization of the conjecture in the case of anisotropic SMA. Consequently, a new equivalent stress and associated strain should be defined. In this context the micromechanical model will be a real help.

  FIG. 1 -Yield Surface in σ 11 ,σ 22 ( ) plane.

FIG. 5 -

 5 FIG. 5-Transformation surface in σ 11 ,σ 22 ( ) plane.


Sébastien Personne

Jérôme Antoni

Jean-Daniel Chazot

Cylindrical cyclic acoustic

come

Introduction a. Objective

Considering the development of a rotating machine, reduction of noise radiation is more and more a topic of concern. This requires the identification and localization of sources of noise.

A useful framework has been developed for signals coming from rotating machine (see [START_REF] Antoni | Cyclostationary modeling of rotating machine vibration signals[END_REF]) that can be modeled as

) ( ) ( ) ( t n t p t x + = , (1)
with ) (t p a deterministic part which can be correlated with a cyclic frame event and ) (t n a random part (for example friction, etc…). The context of concern in this paper is a product development not a failure analysis, thus the periodic part is assumed to be dominant compared to the random part. The goal is to localize the periodic sources radiating noise in space and time (what the author call "cyclic events").

In this context acoustical imaging is an efficient tool to localize sources of vibration from pressure measurement. There is a lot of methods which are unified by the author of [START_REF] Antoni | Focalisation Bayésienne: une approche unifiée du probleme inverse en acoustique[END_REF] with the Bayesian formalism. For a planar geometry, the theory is strongly explained and the efficiency has been well observed [START_REF] Antoni | A Bayesian approach to sound source reconstruction: Optimal basis, regularization and focusing[END_REF]. As far as the authors know this process has not been built for a cylindrical shape; the following article will focus on this point. Readers can find a specific case corresponding to the Statistical Optimized Near-Field Acoustic in cylindrical coordinate in [START_REF] Cho | Source visualization by using statistically optimized near-field acoustical holography in cylindrical coordinates[END_REF]. Usually, specificity of cyclostationary signals is not taken into account and this entails information loss. To overcome this problem, the Authors in [START_REF] Wan | Near field acoustic holography(NAH) for cyclostationary sound field and its application[END_REF] have developed the Near-field acoustic holography for cyclostationary sound field and the source reconstruction is presented in the frequency domain.

The Authors will use another strategy taking advantage of filtering the periodical part of the signal: the Fourier transform can be applied on this part of the signal, apply the image processing in the frequency domain and then use inverse Fourier transform to come back in cyclic time domain. This leads to reconstruction of the cyclic source field in space and cyclic time. Figure 1 illustrates the global procedure. This paragraph is a very short introduction to Cyclostationarity and its properties. For more details on the cyclostationary analysis the reader can refer to [START_REF] Antoni | Cyclostationary modeling of rotating machine vibration signals[END_REF], [START_REF] Antoni | Cyclic spectral analysis in practice[END_REF]. Cyclostationarity gives a theorical framework to processes whose statistical properties are periodically varying in time. By definition it embodies a class of non-stationary stochastic processes, with stationary and deterministic periodic processes as specific cases. When defined on the first order it means that the signal contains some periodic component then CS1 will refer to the deterministic part of the cyclostationary process; on the second order it means the signal may be random while energy is flowing periodically with time.

A formal definition can be done in terms of the expectation operator. A stochastic process, ) (t x , this process is said to be cyclostationary of order i and of period T if any of its i-th order statistics is a periodic function of time with main period T. In this work only the first-order cylostationarity will be considered. Condensing the expression with expected value gives,

{ } ) ( ) ( t x E t m x =
exists and is periodic of period T, i.e.

+ =

(2) In practice the only way to estimate the expected value is to perform a time synchronous averaging with the cycloergodicity assumption.

A stochastic process, ) (t x , is said to be cycloergodic of order 1 and period T if its time-synchronous average is

deterministic, i.e. ) ( ) ( 1 2 1 ) ( lim t kT t x K t x x K K k K T µ = + + = ∑ - = ∞ → . ( 3
)
With the assumption of cycloergodicity, we have

{ } T t x t x E ) ( ) ( = .
(4) Equation ( 3) means that the time synchronous average is an estimator of the expected value. Figure below shows three realizations of a cyclostationary process (left part) and the corresponding time synchronous average (right part), reader can see that the average reveal periodic phenomena. This points out the utility of the cyclostationary framework: since by construction the expected value is periodic, it makes possible to apply in a second step all the powerful tools developed for stationary signals such as the Fourier transform and spectral analysis.

Initially developed for communication signals, this approach seems well adapted for mechanical rotating machines. Due to the kinematics of machines, the resulting acoustic and vibration signals are periodic with respect to some angle of rotation and are therefore cyclostationary with respect to angle. For steady-state operating conditions, angle and time can be reversed and signals are cyclostationar with respect to time as well. The extraction of the CS1 part of the signal will be the first step of the procedure. Those periodic signals will be used for cyclic sound imaging.

c. Bayesian approach to the inverse acoustic problem: optimal basis and regularization

Localization of acoustic sources using discrete measurements of the acoustic field is a recurrent acoustic inverse problem. Many approaches exist, such as Beamforming, Near-field Acoustic Holography, Inverse Boundary Element Method, Helmholtz Equation Least-Squares, etc…Choosing the right technique depends on the acoustic source shape, the array geometry and the frequency range of analysis. However, most of the methods share a process to solve the inverse problem: 1) the source field is reconstructed as linear combination of spatial basis functions, 2) the coefficients of the basis function are computed so as to match the measurement field at the microphones positions. Let s be a radiating source in the acoustic domain and let an array of M microphones measure the acoustic field. The Green function G between the source and the measurement array is supposed to be known analytically or numerically (n refers to the noise of measurement).

n d r G r s p + Γ = ∫ Γ ) ( ) ( . ( 5
)
The aim of the inverse problem is to reconstruct the continuous sound field s(r) from the M discrete measurements. The common formalism to get the reconstructed source ) ( ˆr s can be expressed as:

∑ = = x k k k r c r s 1 ) ( ) ( ˆφ . ( 6
)
where the unknowns of the problem are the basis functions k φ , their coefficients k c and the x dimension of the basis.

The following will give a brief theorical background of the Bayesian approach used to solve for these unknowns -readers should refer to [START_REF] Antoni | Focalisation Bayésienne: une approche unifiée du probleme inverse en acoustique[END_REF] for more details. Bayesian inference is used to solve the problem by considering the unknown quantities as random variables that produce a random source field and by seeking the probability distribution [ ]

p c s \ ) ,
( φ conditioned to the observation of the measurements in the vector p. This is the so-called posterior probability distribution, which may be viewed as a cost function whose maximization will lead to the optimal parameters c ˆand φ ˆ that best explain the measured data.

By using the Baye's rule:

[ ] [ ][ ] [ ] P c s c s p p c s ) , ( ) , ( \ \ ) , ( φ φ φ = , (7)
Where [ ]

The prior distribution [ ]

) , ( φ c s results from any knowledge the experimenter has on the source field before the experiment. A prior can be a spatial information about the radiating regions where the noise is more likely to originate. This information is introduced by an "aperture function " ) ( 2 r S σ that takes positive or zero values on the source surface. This leads to the following structure for the spatial covariance function of the random source field:

{ } ) ' ( ) ( ) ' ( ) ( 2 2 * r r r r s r s S - = Ε δ σ α , (8) normalization 1 ) ( ) ( 2
∑ = - = Ω M k k k k N S U r r G r 1 * 2 / 1 * 2 ) ( ) ( ) ( ψ λ σ , ( 9
)
where G is the vector containing the green functions from the source plane to each microphone of the array, k λ are non-negative scalars, ) (r k ψ eigen-functions over the surface of the source, and k U eigen-vector such that: 0 ...

2 1 ≥ ≥ ≥ ≥ M λ λ λ , ( 10
) kl S l k r dS r r δ σ ψ ψ = - Γ ∫ ) ( ) ( ) ( 2 , ( 11
) kl l k U U δ = * . ( 12
)
It is proved in [START_REF] Cho | Source visualization by using statistically optimized near-field acoustical holography in cylindrical coordinates[END_REF] that the optimal basis functions sought in ( 5) are the eigen-function of (8), i.e. represents the noise-to-signal ratio and plays a crucial role for regularizing the problem.

P U c N k k k k 2 / 1 * 2 2 - Ω + = η λ λ , ( 13
)
To sum up results, the reconstructed source field is ultimately given by (13):

P U r r s N k k M k k k 2 / 1 * 1 2 2 ) ( ) ( ˆ-= Ω + = ∑ ψ η λ λ . ( 14

Cylindrical cyclic acoustic imaging theory a. Scanning measurement for CS1 cyclostationary sound field construction

The shape of the acoustic source is a cylinder of radius r s . The hologram surface will also be a cylinder of radius r h >r s . The goal of the procedure is to measure and reconstruct the cyclic part of the sound field by measuring the pressure on the hologram surface by scanning a virtual array of M microphones, M sequential measurements.

The scanned signals are not temporal synchronous, yet by construction each post processed CS1 signal is phased with a common angular event; therefore the M CS1 signals are synchronous and can be assembled in a "virtual array" of microphones.

b. Bayesian approach to sound source reconstruction on a cylindrical shape

The expression for the radiation of a cylinder is taken from [START_REF] Williams | Fourier Acoustics sound radiation and nearfield acoustical holography[END_REF] ∑ ∫

∞ = -∞ = ∞ ∞ - = m m z z ik im z m dk e e k r P z r p z θ π θ ) , ( 2 1 ) , , ( , (15) where )
z s m s r m r m z m k r P r k H r k H k r P = , ) ( ) ( ) , ( 1 1
where

1 m
H is the mth-order Hankel function and (20) In the computer implementation those infinite summation have to be sampled and truncated.

     < - ≥ - = z z z z r k k for k k k k for k k k , , 2 2 2 2 , (17)
Note that in theory this propagator is calculated for an infinite cylinder and is not valid for a finite cylinder. However, the Baye's formalism allows introduction of a prior in the form of finite support-set aperture function along the z axis that forces a finite cylinder height. The narrower the aperture function and the more localized the prior; how the optimal spatial basis changes accordingly is shown in the following figure. Figure 5 shows that the optimal basis strongly depends on the choice of the aperture function. Man can notice that the spatial basis corresponding to a uniform aperture function (line1 in figure 5) is the specific case use for SONAH in cylindrical coordinates [START_REF] Cho | Source visualization by using statistically optimized near-field acoustical holography in cylindrical coordinates[END_REF].

Numerical simulation a. Source signal simulation

The simulation tries to represent the behavior of a rotating machine; its housing has a cylindrical shape and is the rotating part of the machine. The rotation is due to kinematics events which cause internal forces. Each cyclic event produces a source with its localization and from an external point of view this corresponds to moving sources within the cycle.

Consider a cylinder in motion with a constant rotation speed. Two acoustic sources rotate relatively to this cylinder. For the simulation the CS1 sound field is considered only; accordingly the source variation is made periodic with one rotation but frequency change within this period. The speed is 120rpm and Fs = 2kHz (sampling frequency); From the fix point of view of a sensor on the array, man see two acoustics sources following a circular trajectory as shows the figure 7. Taken the Lagrangian point of view we see the accuracy of the time varying reconstruction of the process. Figure 9 shows the evolution of the particle velocity as if the sensors stuck to each source during motion.

Conclusion

Present work gives a framework for the analysis of a rotating machine with housing as rotating part. This will take advantage of cyclostationary properties of signal generated by rotating machine to localize the acoustical sources in space and time. A simulation demonstrates the ability to reconstruct the sound field over a cylindrical surface from a periodic field measured in cylindrical surface surrounding the source.

This simulation has to be completed with a test on real source to evaluate the global process: from scanning measurement to sound field reconstruction. Particularly the influence of the phase error due to the angular resampling process, see [START_REF] Bonnardot | comparaison entre les analyses angulaire et temporelle des signaux vibratoires de machines tournantes, etude du concept de cyclostationnarité floue[END_REF]. Reader can see in [START_REF] Boucaud | Improved acoustical sources location in gear power transmissions[END_REF] that, measurement done by step and synchronized by a precise angular event enable to find an accurate source location in a gear power transmission. Author's work will be to follow the same way with a "Cylindrical robot" for a cylindrical source.

FIGURE 1 :

 1 FIGURE 1: synopsis of the procedure b. Cyclostationary: definition and properties

FIGURE 2 :

 2 FIGURE 2: Typology of stochastic processes

FIGURE 3 :

 3 FIGURE 3: illustration of the time synchronous average and expected value construction

2 β

 2 function, reflecting the probability of observing p given the source field priori probability distribution of the source before the data are measured, and [ ] P is the evidence.It is important to notice that with this formulation, the inverse problem (inferring s(c, Φ) from p) has been turned in its direct problem (inferring p from s(c, Φ)). probability distribution as the measurement noise n (n is the vector of noise measurement for the array). According to the central limit theorem applied to the Fourier transform, the measurement noise in frequency domain is a circular complex Gaussian random variable. Introducing the covariance matrix reflects the mean energy of the noise.

2 α

 2 reflects the mean energy of the source field.With the expression of the probability distributions, the inverse problem corresponds to the maximization of the main results of the algorithm are reported. We define the singular value decomposition of the propagation operator G through the aperture function and whitened by the covariance matrix of the noise N Ω as follows:

  . The optimal coefficient k c sought in (5) are found as

  )In the present work, the reconstructed source is the CS1 part of the sound field. So the vector p in (13) contains the CS1 part of the M signals measured by the microphones of the array.

FIGURE 4 :

 4 FIGURE 4: continuous source surface and discrete position for the microphones of the virtual array

  θ , the m-th circumferential component of the sound field. The wavenumber spectrum at radius r can also be expressed in terms of the wave number spectrum of the sound field on a cylindrical source of radius r s :) , (

with k= c ω ,

 ω ω the angular frequency, c the celerity and the convention for time dependence factor is

σ

  Module of the 10 first eigen functions of the propagation operator G through the aperture function for f=500HzPrior is only in the z direction: uniform localization of sources around the cylinderWindows is normalized as follows Without prior information, i.e. a uniform window for the aperture functionWith a prior information given by the aperture function With a prior information given by the aperture function

FIGURE 5 :

 5 FIGURE 5: Evolution of the basis function for source reconstruction with prior on the source location along Z axis.

FIGURE 6 :

 6 FIGURE 6: Evolution of sources with respect to time and associated instantaneous frequency. Signals are shown over 1 cycle, T cycle =0,5s.

FIGURE 7 :

 7 FIGURE 7: geometry used in simulations b.Source reconstruction by Bayesian approachThe Eulerian point of view is used to obtain a clearer idea of the source localization at one particular instant. Gathering successive pictures give a movie of the source motions. Figure8presents two pictures taken for sample 100 and sample 900.

FIGURE 8 :

 8 FIGURE 8: comparison of simulated source field and the reconstructed source field for two specific time in the cycle

FIGURE 9 :

 9 FIGURE 9: reconstruction of each in time