Nicolas Simond

Patrick Rives

What can be done with an embedded stereo-rig in urban environments ?

Keywords: stereo-vision, mobile robotics, urban environment, super-homography, trajectography, 2.5D reconstruction .2.10.4: Motion I.4.6.0: Edge and feature detection

Introduction

A prerequisite to the design of future Advanced Driver Assistance Systems for cars is a sensing system providing all the information required for high-level driving assistance tasks. In particular, such a sensing system should be able to provide a robust and accurate real time localization in constrained urban environments. Moreover, if we want to address Autonomous Guided Vehicles (AGVs) issues, we must also take into account on board sensors able to provide representations of the environment supporting obstacle-free trajectory planning algorithms and controller design based on reactive sensor-based control tasks.

Commonly, the localization in outdoor environment is based on the hybridization of GPS with proprioceptive sensors embedded on the vehicle like odometers and/or Inertial Navigation System (INS). In urban environment, this process is enhanced using a prior knowledge from maps thanks to a map matching algorithm which constraints the trajectory of the vehicle to belong to the road. Nevertheless, the success of such techniques is closely related to the reliability of the GPS data and its availability during the navigation task. Practically, in a dense urban environment (ie. old city downtown), these prerequisites hold rarely. Due to the presence of vertical structures in the architecture, the visibility of a sufficient number of satellites for positioning (at least 4) is time dependent and the signals are frequently corrupted by multipaths propagation. When reliable GPS data are missing, the localization process is only supported by a dead reckoning method based on the integration of odometry data provided by internal proprioceptive sensors. It is well-known that such sensors are subject to an important drift leading to a biased localization.

On the other hand, it is now clearly admitted that, in a near future, every car will be equipped with on-board exteroceptive sensors like vision, lidar or radar for ADAS applications. Such sensors are well-adapted to capture the interactions between the vehicle and the local environment in the vehicle-centered frame. Under certain assumptions, from these data, it is possible to simultaneously estimate the local structure of the environment and the ego-motion of the vehicle. This problem often called Simultaneous Localization and Mapping (SLAM) problem was intensively studied in the case of indoors robotics. Nevertheless, the transposition of SLAM indoors methods to applications in urban environments is not straightforward due to the complexity and the dynamical character of the urban context. SLAM algorithms for specific inner-city applications have been developed for the last five years. To ego-locate from the cluttered environment, the mobile robots have first to segment out the moving objects before building maps and localization task. Hence, the CMU Navlab purposed in 2003 to accomplish SLAM with Detection And Tracking of Moving Objects (DATMO) [?] that became in a unified approach SLAM and Moving Object Tracking (SLAMMOT) [?], in 2004. When a model of the environment is known, the vehicle localization tends to a Structure From Motion (SFM) problem ([?]) because opposite to SLAM methods, the duration of the environment visibility is very short with a traveling vehicle. Recently, collaborative GPS/INS navigation ([?]) is investigated to provide data localization to a vehicle in a GPS shading area.

Basically, on-board vision was implemented to improve the robustness of the navigation task by analyzing the evolution of the environment. Vision has in fact many applications on obstacle detection and lateral control of a path following limited by boundaries (line or lane tracking). Nevertheless, many authors try to estimate the ego-motion of the camera(s) assuming a prior knowledge on the environment:

• optic flow ([?,?]) is used in case of unstructured environment or textureless environment, • direct method ([?,?]) provides reliable results assuming a calibrated camera and a motion model reduced to the three main parameters, • homography ([?]) to detect planes in the environment and compute then the relative motion of the camera.

Note that the three last references allow the estimation of a vehicle ego-motion following a textureless road.

Due to the fact that urban environment are structured, the localization of a vehicle is affordable assuming a transposition of visual methods dedicated to indoor environment. The relative camera orientation can be computed from the location of vanishing points ([?,?]). On the opposite way, the navigation of a mobile robot requires the fusion of data from odometer, INS, GPS, vision and an a priori knowledge of the static environment as a Geographical Information Systems (GIS) to improve the reliability of the global localization with identifying "natural" beacons like [?] of the NYU Avenue Project. The introduction of artificial beacons in a city is thus not an acceptable solution due to the scale problem. Authors prefer identifying vertical edges which are precisely referenced in GIS ([?,?]). Another way consists on the update of the current view with a set of images, recorded and geo-referenced during a calibration phase ([?,?]).

We present in this article a novel method to estimate the trajectory of a vehicle navigating in an urban environment thanks to the data provided by an onboard stereo-vision system. We assume a canyon-like environment which contains some static planes (road, frontages) where characteristic features (points, lines) can be extracted. The purposed method copes with the dense traffic conditions: the free space required (first ten meters in front of the vehicle) is slightly equivalent to the security distance between two vehicles. Moreover, the estimation of the pose of the stereo rig along the navigation path allows us to build a polyhedral model of the environment using a SLAM technique.

The paper is organized as follows. In the first section, a simplified model of the urban environment is proposed and we present an overview of the method.

The second section focuses on the segmentation of the main planes along the sequence of images. The third section is devoted to the robust estimation of the coplanar features which belong to the road plane from the sequences of images provided by the stereo rig. In the fourth section, we develop a method based on the decomposition of the homographies linking the different views in order to estimate the ego-motion and the geometry of the planes in the scene. Finally, thanks to a mosaicing-like technique, we compute both the trajectory of the vehicle and a representation of the underlying road. Experimental issues from real data are presented and discussed. We conclude this article in the last section by proposing perspectives.

Problem Statement

Viewed from inside the vehicle, urbanized environment can be viewed as a kind of canyon where the ground is bounded by the vertical frontages of the buildings. The scene can also be modeled by a hallway model where each main surface is assumed to be locally planar at the foreground. The frontages are generally formed with sets of quasi vertical planes containing patterns whose edges are generally aligned with two main directions that are the road and the vertical axis. Concerning the road, signs lying on the ground, they are mainly aligned with the longitudinal and the transversal road axis. We consider the vehicle is equipped with an on-board stereo-rig providing standard gray level images at video rate. The vehicle does not exceed the speed limit in urban environment (abroad 50 km/h in European countries). Consequently, we also assume small displacements of the features between two successive images in the sequence. We denote the discrete time index by the variable k. At each iteration (k), two [n raw × n col ] images I k l and I k r (respectively, left and right), are recorded.

Mathematical prerequisites

Let P s be a 3D point of the environment, its projection onto the image is represented by a [3 × 1] homogeneous vector p k s = [u, v, 1] t . The contour L s in the 3D environment projects onto the image as a [3 × 1] homogeneous vector l k s = [cos(θ), sin(θ), -ρ] t where θ and ρ are the polar representation of the line in the image frame.

A 3D point in the scene is related to its projection in the image by a projective equality: p s ∝ K[R,t].P s where K represents the [3 × 3] matrix of intrinsic parameters, R and t the rotation matrix and the translation vector between the scene and the camera frames. When a subset of features belong to the same plane, their projections in image plane are also related to the linear combinations between their coordinates in the original plane.

We hence consider different points of view of a same plane like in Fig. 1. It is now wellknown [?,?] that for each couple of images I a and I b , it exists a mapping from a projective plane to a projective plane, called homography, defined up to a scale factor (3 × 3 matrix but only 8 independent parameters). This homography links the coordinates p a and p b of the features in the images I a and I b such that p b ∝ H ba p a . According that each corresponding points (resp. line) in a couple of images resulting from the projection of a unique point (resp. line) in the scene provides two independent equations of constraint, only four independent points (resp. lines) are required to compute an homography. In the Fig. 1, homographies induced by the π plane between the stereo images are noted H st , whereas homographies between consecutive left images or right are noted H lr .

The estimation of the camera motion (rotation and translation) [R ba , t ba ] between the images I a and I b from the homography H ba , requires the knowledge of the calibration matrix K a , K b , the unit vector normal n a and the distance to the plane d a , both expressed in the first camera frame:

H ba = K b R ba - t ba .n t a d a K -1 a (1)

Overview of the Method

Among the different planes constituting the urban canyon, the road plane takes a central importance. Sometimes, in large avenues or in presence of trees, it will be the unique plane which remains visible when the vehicle is moving. So, it is fundamental for navigation and guidance purposes to robustly and accurately detect and characterize the road plane. Unfortunately, due to the uniformity of the surface, the lack of texture and the few of characteristic features (excepted the road signs), this detection will be the most difficult.

In the next two sections, we detail how we can estimate the ego-motion of the vehicle using the features extracted from the road plane only. The results will be extended to the others planes when they are available, in the fourth section.

The basic idea of the method is to take advantage both the spatial constraints linking the left and right cameras to extract the features and the temporal constraints between the successive views. An overview of the complete algorithm is depicted in the Fig. 2. The estimation of the stereo-rig egomotion can be splitted in four linked sequential stages:

(

Segmentation of the road plane

In an urban context, the road is commonly viewed by the on-board cameras as the dark homogeneous region limited by the light road markers, located on the bottom of the image. Assuming that some of the road markers are visible in the foreground, they can be used to estimate the road direction in order to segment the road area in the images. The identification of the road markers remains nevertheless complex due to presence of kerb-sides, intersection areas or the obstruction by dynamical obstacles. Consequently, using only a data-based approach for segmenting the road is not robust enough and it will be benefit to introduce some global constraints from a model. Let us assume i) the sides of the road are parallel and can be locally approximated by segments of lines at the foreground, ii) most of the edges in a canyon-like urban environment are aligned with three main directions: the longitudinal and transversal axis of the road and the vertical axis.

Extraction of the Dominant Vanishing Point

Under a perspective projection, a set of parallel lines in the 3D scene projects onto the image as a pencil of 2D lines which converge to the same point, called vanishing point. The vanishing point is characteristic of a direction in the scene and it belongs to the plane at the infinity. As the projective homography linking the plane at the infinity and the image plane depends only on the rotations, the vanishing point will be invariant with regard to a translation motion of the camera. Such a property will be useful to detect and track the vanishing point corresponding to the direction of the road axis referred below as the Dominant Vanishing Point (DVP) x.

Considering our embedded stereo-rig, the variation of the DVP coordinates has two possible origins: the variation of the pan angle when the vehicle is turning and the variation of the tilt angle when the vehicle is accelerating or braking.

Assuming the baseline of the stereo-rig is parallel to the road plane and the tilt angle remains small, the model of evolution of the DVP coordinates can be selected as follows: a constant velocity model for abscissas and a stationary model for ordinate.

Thanks to a Kalman filtering technique, the estimation and the tracking of the DVP during the vehicle motion can be performed robustly. The observation vector Y is computed from a subset of lines S = {l t 0 . . . l t M } built from the segments extracted into the current image using a Canny operator. To be candidate the lines have to converge in a neighborough of the predicted DVP. The observation vector Y is then computed as the image point which minimizes the weighted sum of distances to the subset of selected lines.

Y = min x M s=1 (w s .l t s .x) 2 (2)
with w s a weight that represents the sum of the length of the segments contributing to the line l t s .

To be exhaustive on this topic (see [?] for more details), we would like to stress that the quality of the estimation of the DVP is closely related to the solvency of the equation ( 2). When all the lines have quite similar orientations in the image (like in Fig. 3), the solution of the equation ( 2) becomes very instable. In this case, to avoid discontinuities in the DVP estimation due to a lack of constraints, we propagate the prediction of the ordinate provided by the Kalman filter instead of the estimated value. The Fig. 4 illustrates the result of the DVP detection and tracking algorithm applied to a stereo video sequence involving an overtaking maneuver [START_REF] Antone | Automatic recovery of relative camera rotations for urban scenes[END_REF] . The vehicle follows a straight path, stops before a parked vehicle, overtakes it then continue its straight motion.

The plots represent the evolutions of the DVP coordinates estimated by the Kalman filter in the left and right images. The smoothness of the chronograms confirm us on the adequation of the models. A DVP abscissa constant corresponds to piece of straight path. The ordinate variations due to the pitching appear limited in comparison of the abscissas ones that justifies the stationary model of evolution. The small variations between the frames [960;990] results from the braking and acceleration phasis during the overtaking maneuver.

Characterization of the traffic lanes.

Most of authors assume that the road limits can be easily detected due to their artificial (road markers) or natural (kerbs, soft shoulders or green strips, sidewalks) contrast with the pavement. The detection of these road markers is intensively used in many ADAS applications or for lateral control in automatic guidance schemes [?,?,?]. Most of the time, urban roads are structured: lanes of navigation have constant width and are highlighted with painted lines. But in standard traffic conditions, only the lane where the vehicle is moving on, can be expected to be observed by the cameras, the others lanes are almost always occluded by the flow of vehicles and can be perceived from time to time only.

¿From these reasons, we state the road detection problem as a parameter estimation problem. We assume a region of interest corresponding to the road in the image is defined by the current DVP coordinates and two points for a one lane road (resp. three points for a two lanes road) resulting from the intersection of the median axis of the left and right road markers for a one lane road (resp. the left, middle and right road markers for a two lanes road) with the bottom of the image, like in Fig. 5.

At each abscissa of these points is associated a constant velocity filter and the abscissa difference between two medians abscissa is assumed to be constant. The observation of the filter is computed from the detection in the image of stripes which satisfy a dark-light-dark gradient transition and locally oriented towards the DVP. Moreover, more importance are accorded to the stripes detected close to the bottom of the image (just ahead the vehicle). At the initialization step, a priori knowledge on the number of lanes and the the painted lane width (in teens of pixels) is provided. In the first couple of images, the stripes have to be pointed out manually by a user if the automatic tracking fails.

During the tracking process when a stripe jumps outside the limits of the image thanks to the constant lane width constraint, it remains possible to perform the abscissas prediction of the loss stripe median. That is the case for two of three stripes (middle and right) between frames [936;952] in the Fig. 6, when the test vehicle has large rotation motions to overtake the parked car. The quality of the estimated abscissas allows the identification of the two others road markers as soon as they go back inside the field of view.

Road plane Estimation

The segmentation stage described above, can be viewed as a pre-processing stage which aims at drastically reducing the search area for the features liable to belong to the road plane. Conversely, due to the presence of dynamic obstacles, all the features inside the Region Of Interest (ROI) are not really lying on the road plane. In this section, we present a novel method based on the computation of the homographies induced by the road plane in the different views, for classifying the features extracted from the ROI as features lying or not on the road plane.

In our particular application, a first idea will be to use the DVP point and the road markers previously extracted for computing the homography. Unfortunately, such a choice is inadequate due to the fact that the projection of parallel lines in the scene yields to linear dependent lines in the images. Others features are needed and we propose to exploit also the points of interest extracted from the ROI.

Computing the Homography

Two types of homography are computed at each iteration k: the stereo homography H k st between the left and right images of the stereo-rig and the homographies H k l (resp. H k r ) between two consecutive left (resp. right) images. The stereo homography allows the extraction of the coplanar features into the ROI assuming that a reliable prediction is avalaible while the computation of H k lr , which are our goals, are performed when the majority of the coplanar features are known in both images.

The road pavement appears as a textureless region where the use of SIFT keypoints ([?]) is not well-adapted: the matching of keypoints extracted in this area between two stereo images provides unreliable results. Therefore we use an Harris detector which extract most of the stable feature points (FPs) from the road signs (lines, stripes, arrows, letters, etc.) or from the static and dynamic obstacles present on the road. Unfortunately, only a few subset of FPs extracted from the obstacles do represent an intersection with the road plane: the others should be discarded at the processing level.

Homography computed from stereo-vision

Finding correspondences between FPs: The purposed method is based on a variant of the Iterative Closest Point (ICP) algorithm [?] that is an alignment algorithm that works in three phases : 1) establish correspondence between pairs of features in the two structures that are to be aligned on proximity, 2) estimate the rigid transformation that best maps the first member onto the second and 3) apply transformation to all features in the first structure. The main difficulty we face is the matching of the two sets of coplanar features which are not sufficiently characterized and uniformly distributed to compute the induced homography.

The Scott Longuet-Higgins ([?]) method consists on the SVD decomposition of a proximity matrix Dist which contains the inter-distance between the two clouds of FPs, weighted by a gaussian.

Dist(i, j) = exp - dist 2 (p i ,p j ) 2.σ 2 (3)
The SVD decomposition verifies Dist = U.D.V t , the replacement of the diagonal elements of D by an identity matrix I provides a new proximity matrix Dist svd = U.I.V t which both respects the exclusion and proximity principles. Each element of Dist svd belongs to the interval [0;1]. Couples of FPs which maximizes simultaneously raw and column of this particular matrix verify the one to one exclusion principle according to the distance criterion.

The algorithm of Pilu ([?]

) is dedicated to the matching of stereo images (small baseline). It performs the last method with purposing an appropriate correspondence strength matrix which integrates both a geometric and an intensity relationship. Parallel to the Dist matrix, the Corr matrix contains the correlation between gray level values of patchs centered on each FP. Due to elements of Corr are normalized between [-1;1], the correspondence strength matrix is:

L(i, j) = 1 + Corr(i, j) 2 .Dist(i, j) (4)
All the couples which maximize L svd, SVD decomposition of L, are not correctly matched. The reliability of the method depends on the threshold we fix to discard couples which can be considered as outliers.

A priori knowledge on the cameras motion: Assuming a rigidity constraint between the two cameras of the stereo-rig, the identification of the road plane and the extraction of coplanar features become easier. The rotation and trans-lation motions (R st , t st ) between the left and right cameras are then fixed and the stereo homography induced by the road plane only depends on the parameters of the plane (n k st , d k st ) expressed in the left (right) camera frame. Furthermore, the high frame rate of the video sequence allows us to consider small variations in the parameters of the road plane between two successive image acquisitions.

n k st ≃ n k-1 st and d k st ≃ d k-1 st (5)
Consequently, reliable predictions of the stereo homography can be done between two iterations using a simple stationary model:

Ĥk st = H k-1 st (6)
Using this prediction, we can search for potential matches between the stereo couple of images. In practice, the hypothesis of a stationary model does not hold due to the perturbations. Thanks to a recursive estimation algorithm, we refine a posteriori the estimation of the stereo homography which allows us to detect and discard the wrong matches corresponding to FPs which are not really lying on the road plane.

Adaptation of the Pilu method to our case: Assuming the homography prediction, we also search for the correspondences between the two sets of FPs extracted from the left and right current images. A couple of FPs is potentially candidate to the match if and only if it satisfies the crossed re-projection criterion, like in Fig. 8 :

     dist(H k st p ls , p r j ) < th dist dist(p ls , (H k st ) -1 p r j ) < th dist (7)
All couples of candidates which do not satisfy the distance criterion, are discarded. The second selection rule consists on computing the normalized correlation between gray level of regions [21 × 21] pixels centered on each couple of FPs candidates. Only the couple of FPs satisfying a score greater than th corr = 0.5 (empirically fixed) are finally labeled as candidates for the matching process. The strength matrix is also sparse that speeds up the SVD computation and reduces the probability of outliers.

Refining the homography estimation: The quality of the stereo homography computation depends on the spatial distribution of the true coplanar features and on the capability to discard the FPs which are detected on the obstacles over the road plane. The stereo homography computation can hence be resumed as a recursive estimation process which toggles: (1) the selection of potentially coplanar FPs assuming the current estimation of the homography as described above, (2) the update of the homography using the couples of FPs matched successfully,

The process stops when the list of selected couples remains stable between two iterations. That is generally the case after two or three iterations ( We would like to stress that the computation of the correlation score between squared patchs centered on FPs lying on a plane which is nearly parallel to the camera axis is not optimal due to the perspective distortion. Assuming that a homography prediction is available, a solution consists on computing the warping of the ROI with the homography prediction and then compute the correlation score. The experiments show that when the FPs have not a uniform spatial distribution, the recursive estimation of the homography sometimes diverges.

Initialization step: At the beginning of the sequence, the prediction of the stereo homography is unknown. Assuming that the majority of the detected FPs in the ROI belong to the road plane, we compute a first estimate of the homography using a RANSAC method ([?]). The stereo homography is computed from a set of four couples of FPs (verifying the correlation criterion) selected randomly. The number of random selections is determined assuming a fixed rate of outliers in the sample.

Homography computed from dynamic vision

The stereo-vision approach provides us with a snapshot where all parts of the environment look like static. The homography computed from stereo-vision gives us an information of the geometry of the environment (i.e. partitioning features lying on the same planes). Complementary, the left and right video sequences informs us about the dynamic parts of the environment and about the motion of the vehicle. However, only just observing the optical flow field between two successive images, it is not obvious to separate the component due to the camera ego-motion (i.e. induced by static objects) and the component relative to the dynamic objects in the scene. To overpass this difficulty, we propose to track in the video sequence only the FPs satisfying the stereo homography of the road plane.

With regard to the computation of the stereo homography, the difficulty we have to face is the lack on knowledge about the vehicle motion. Consequently, we can not compute a prediction of the homography linking two successive images. As done by several authors [?,?], we assume a forward motion of the vehicle composed by a translation along the longitudinal axis, a rotation around the vertical axis (pan angle) and a rotation around the lateral axis (tilt angle).

These approximations allow hence the selection of corresponding candidates for each FP selected with the stereo homography in both images. Couples of matched FPs are also extracted with the advanced Pilu algorithm. The homography is then computed with the same method as the stereo homography computation. At this step, the conditions are identical: we benefit from couples of matched FPs, a prediction of homography and characteristics of the road markers. We present in Fig. 10 some results of homography computations in basic cases.

Integration of spatio-temporal constraints

The stability of the homography computation depends on the coplanar features we achieve to detect. Whatever the type of computed homography may be, some couples of FPs located over the road plane are erroneously selected. We aim at rejecting these outliers for improving the quality of the homographies. The basic idea is to exploit all the constraints linking a larger set of views assuming that the composition rule H ca = H cb H ba is verified whatever the images I a , I b and I c .

The concept of super-homography

Malis and Cipolla [?] describe an efficient method to impose the constraints between the homographies computed from a sequence of views of a planar structure. Tacking into account multiple views provides a set of constraints between the coordinates of coplanar features in different views. Moreover, it minimizes the effects of errors on the matching step and reduces numerical instability when the motion between two views is not enough significant.

Such a method improves the consistency of the current homographies (H k st , H k l , H k r ) estimation. Nevertheless, we also face a hard compromise between increasing the distance between two views to reduce the numerical instability and keeping a significant number of matched features to constraint the homographies computations. The method introduces a super-homography matrix H k [3m × 3m] where m is the number of views. This new entity contains all the [3 × 3] homographies between the (m(m -1)) couples of different views. We personally use three stereo couples corresponding to the frames (k -2), (k -1) and (k), the number of images is then m = 6. Obviously, solving the constraints imposed by the super homography requires the observation and the tracking of the FPs in the (k -2), (k -1) and (k) couples of images. However, as we show below, the presence of an FP in the whole set of m images is not necessary and only a partial observation in a subset of m will be sufficient in practice.

The super feature points

At the preceeding stages, the computations of current homographies (H k st , H k l , H k r ) provide three lists of correspondences between features of the last two couples of stereo-images. The identification and verification of the cross-links between the couples of FPs is processed with fulfilling a table where each column represents a 3D point and the raws the images of the last and current couples. This table allows us to verify the relevance of the current matching steps between the last four images:

p k rs = [H k st .H k l ].p k-1 ls = [H k r .H k-1 st ].p k-1 ls (8)
Let us now introduce a new entity: the super feature point (SFP) which is a [3m × 1] vector which contains the homogeneous coordinates in the m images of a point lying on a road plane. Generally, several points are not detected all along the views. In this case, the coordinates of unknown projections will be initialized with zeros and estimated during the computation of the super homography.

The virtual feature points

As the previous examples shown, the computation of the homography is highly sensitive with respect to the spatial distribution of the 3D points in the plane.

The pencil of road marker medians is certainly the most reliable feature we succeed to track over the video-sequence. Its detection is robust with regard to the local occlusions induced by the presence of obstacles on the road. Thanks to the property, we can define a new type of features, called virtual feature points (VFPs), to constraint the road plane projection in areas where no SFP is detected.

The VFPs are built from the intersections of the pencil of road marker medians with virtual lines lying on the road, defined by couples of coplanar SFPs observed in the m views. Best results are obtained with virtual lines located in the bottom part of the ROI with different orientations. We personally look for nine lines whose orientations are approximatively distribute between [-45 • ; 45 • ] with a step of 10 • , like in Fig. 11.

The computation of super-homography

The concatenation of the SFPs and VFPs coordinates in the m images composes a global feature matrix 0 F k [3m × n F ] where n F = n F P + n V F P :

0 F k = [ 0 P k 1 , 0 P k 2 , ..., 0 P k n FP , 0 VP k 1 , 0 VP k 2 , ..., 0 VP k n VFP ] (9)
The homographies which form 0 H k are initialized from the composition of homographies estimated at the last and current frames, e.g.:

0 H k (I k l I k-2 r ) ∝ H k l .H k-1 l .H k-2 st ( 10
)
At this point, SFPs and VFPs are both vectors of projections coordinates observed in the m images. The first estimation of each column of1 F k is also computed according to some projection coordinates are unknown:

1 m * 0 H k . 0 P k s ( 11
)
where m * ≤ m represents the number of known projections of the s th SFP with s ∈ [1; n F ].

When the composition rule is verifyied for all homographies, rank(H k ) = 3 whatever m ≥ 3. Due to all the numerical errors, rank( 0 H k ) > 3. Therefore an iterative estimation is required which toggles until q ≥ 0 with the estimation of SFPs coordinates according to the last super-homography computation and the computation of all the homographies induced by SFPs projections in all images:

     q F k ∝ 1 m q H k-1 . (q-1) F k (q-1) F k (M r , :) ∝ q H k (M r , M c ). (q-1) F k (M c , :) (12)
where {M r , M c } respectively represent numbers of raws and numbers of columns relative to the images I mr and I mc for the super-homography with {m r , mc} ∈ [1; m] and m r = m c .

The recursive process generally stops after q = 2 iterations unless the system was ill-conditioned due to errors on feature matching or lack of constraints with the features spatial distribution. The ideal sub-pixellic homogeneous coordinates of the features in the current stereo images can be extracted from q F k while the current homographies can be extracted from the final estimation of q H k :

H k =                  I 3 (H k-2 st ) -1 • • • • • • • • • • • • H k-2 st I 3 • • • • • • • • • • • • • • • • • • I 3 • • • • • • • • • . . . . . . . . . I 3 . . . . . . • • • • • • H k l • • • I 3 (H k st ) -1 • • • • • • • • • H k r H k st I 3                  with P k s = p k-
Note that the upper matrix of H k is the inverse transpose of the lower part that makes the 0 H k computation easiest.

The Fig. 12 shows the utility of the VFPs computation when the FPs spatial distribution is not uniform in the ROI or when obstacles obstruct the ROI clear view. Initially, four of the five matched FPs were aligned with the road marker, the induced linear system was ill-conditioned. The introduction of the SFPs, light '⋄' and dark '⋆' respectively in current and past image, drastically improves the quality of estimation: pencil of VLs correctly superimpose and new couples of FPs are identified. Among them two are not aligned with the left road marker.

6 Experimental Result

Trajectography of the vehicle

Assuming the intrinsic parameters of the cameras are known (cf. 1), the introduction of the super-homography and the VFPs allow the vehicle trajectory calculation along some hundreds of meters. The integration of the homographies H k lr , extracted from the H k over the sequence, provides the current poses of the vehicle in an absolute framework. To retrieve the vehicle motion close to the reality along the Versailles sequence, the original orientation of the right camera is estimated to [2; -6.5; 0] • in the Fig. 13. The sequence represents a forward motion of the vehicle along more than 250 m with a change of lane at the end.

Reconstruction of the road plane

As a straightforward result of the homography-based method, we are able to reconstruct a 2.5D model of the road plane in the scene and warp on it the images taken during the motion (see Fig. 14). The bird eye view of the road made of the warping of 170 images of the ROI segmented in the right images of the stereo-rig. The discontinuity of the road at the middle of the image is due to a lack of texture in the preceding road area. The extracted SFPs do not constraint sufficiently the super-homography which is finally incorrect.

We expect that such a representation could have great uses for the urban policy of the cities: the pavement service, cadastre, real-time navigation systems providers. Direct measures are indeed affordable in hybrid image, that offers new solutions for these users. Each frame is also represented in the first view of the sequence with considering the composition rule of homographies between the first and the current frame.

2.5D modeling of the environment

In the specific context of urban environment, a 2.5D modeling can be performed according to the frontages of buildings are mainly composed with vertical planes. The frontages are structured areas where several corners and edges can be extracted. The method developed on the road plane can hence be easily transposed to the vertical planes. Opposite to the FPs lying on the road plane, the FPs extracted on the frontages have a radial motion with a dominant horizontal motion towards the vertical boundaries of images.

We show in Fig. 15 the result of super-homography computation applied on the three main planes of the scene. The vertical planes are first segmented according to the hallway model applied on the environment after the detection of the road area. The verification of the SFPs planarity allows the identification of new regions in images which are not projections of the segmented planes: they are assumed as obstacles. The obstacle detection is the required task to perform the robust navigation of a vehicle. In a piecewise planar scene, this method provides rapidly reliable results. The localization of the vehicle assuming the super-homography computation from the frontages can not be performed with this configuration of cameras.

The clear view of the frontages is limited by the obstacles and the distance between frontages and cameras reduce their relative motions in images. In other words, the common frontages areas are smaller than the road one whatever the couple of images. The homographies induced by frontages are also less accurate than the homography induced by the road, they could still be used as predictions on the vehicle motion when the road plane is occluded.

Conclusion and perspectives

We detail in this article what sort of informations can be extracted from an on-board stereo-rig. Assuming that urban environments can be assumed as a piecewise planar environment, we present promising results which deals with the autonomous guidance of a vehicle in the urban conditions. The main goal of our work is a vision-based localization method which estimate the vehicle motion according to the tracking of the static scene. Highly urbanized environments can be generally modeled with a hallway model and have structured roads where the lanes are distinguished with road markers. In dense traffic conditions, the clear view of the static scene is limited. We nevertheless develop a method which segment the road area in images to extract coplanar features (lines and points) according to the stereo-vision constraints. The vehicle motion is then estimated with computing the homography induced by features lying on the road plane between two frames.

Dynamic and static obstacles reduce the clear view of the structured static scene. They provide furthermore features which corrupt the matching stage of coplanar features. The homographies computation are also improved by verifying the features coplanarity in multiple views. The super-homography allows it in a single computation and improve the reliability of the coplanar features extraction. The method can be easily transposed to each plane which structured the static scene.

Two material drawbacks limit for the moment the method. First, the field of view of the road is too restricted when the vehicle turns. The use of wide field cameras or the introduction of a third camera certainly allows a widening of the field of view that can maintain an optimal view of the road area in images. Second, the major difficulty we face is to do not access to a prediction of the vehicle motion between two consecutive frames. We plan to synchronize the frames acquisition with an odometer pulse that has the great advantages to allow a reduction of the frame-rate when the vehicle is not running and otherwise to bundle the relative motion of the static environment between images.

In case of post-processing applications, an accurate segmentation of the road plane can be performed assuming homography computations with direct methods ([?]). Direct methods are indeed well-adapted to textureless area but requires rectified images of the camera and a calibration stage which can be provide by our method. Otherwise, the homography parametrization developed by [?] should bundle the variations of the homography parameters and avoid discontinuities in estimated motions. Auto-calibration of the stereo-rig can be performed assuming that the parallelism and orthogonality of some planes which structure the urban environments.

[?] exploit the interest to perform an on-board camera calibration when stairs or pedestrian crossing areas are observed.

Fig. 1 .

 1 Fig. 1. Motion of a stereo-rig in front of a plane between the discrete time (k-1) and (k). Rotation (R) and translation t motions of cameras are constraint by the rigid transformation between the cameras framework.

Fig. 2 .

 2 Fig. 2. Overview of the method.

  1) the segmentation of the road plane in the left I l and right I r images by taking account the Road Limits and parallel lines of the environment, thanks to the extraction of the Dominant Vanishing Point x, (2) the extraction of candidate coplanar Feature Points p in the segmented region which satisfy the spatial constraints induced by the stereo-rig geometry (stereo homography H st ), (3) the candidate coplanar features are tracked in the sequence and the false matches corresponding to moving parts in the scene are rejected, (4) the refinement of the feature selection by using homography constraints linking multiple views in a bounded horizon time: computation of the super-homography H .

Fig. 3 .

 3 Fig.3. Case where the estimation process is not constrained sufficiently.

Fig. 4 .

 4 Fig. 4. Chronogram of the left and right DVP coordinates (abscissas at the top, ordinates at the bottom) along the Antibes sequence.

Fig. 5 .

 5 Fig. 5. Tracked road markers in both images of the stereo-rig in case of a two-lanes road. The dashed lines are others edge-lines which converge to the DVP location.

Fig. 6 .

 6 Fig. 6. Chronogram of the road markers characteristics along the sequence: abscissa and orientation of the medians in both images then the relative width of the stripes measured at the image bottom.

Fig. 7 .

 7 Fig.7. The grid applied into the ROI. Among all the detected FPs '.', only the '+' are selected to compute the current homographies.The accuracy and the stability in the computation of the homography is closely related to the spatial distribution of the FPs into the ROI. The road markers generally provide the majority of the FPs. To improve the spatial distribution, we only keep the most representative FPs in each subregion of a grid applied on the ROI, as in Fig 7.The selection is performed with the Harris score. The road is hence divided into two types of region: those where are located on the road markers and the others. The horizontals split also the ROI with equal height slices that artificially balances the number of FPs detected at the foreground in comparison of FPs detected far from the camera. Some grid areas still remain empty because of the detected FPs have score lower than a critical threshold.

Fig. 8 .

 8 Fig. 8. Locations of the detected FPs ('+') and their projections ('.') in the other image, computed with the homography prediction Ĥk st .

Fig. 9 .

 9 Fig. 9. Warping of I l into ∇ l with H st . On the left, the spatial distribution of the selected FPs is uniform, opposite to the case on the right. Homographies are computing with couples of coplanar features: '+' and dashed lines of I l are in correspondence with 'o' and solid lines in I r .

Fig. 10 .

 10 Fig. 10. Relative motions of the coplanar FPs ('+') according to a straight translation (left) and a dominant rotation motion of the vehicle (right).

Fig. 11 .

 11 Fig. 11. Virtual feature points ('⋄') are intersections of pencil of road markers medians (solid) and the virtual lines (dashed), defined with couples of SFPs, whose projections are known in the m views.

Fig. 12 .

 12 Fig. 12. Influence of the feature spatial distribution to compute homographies. The hybrid images are composed with current image (dark features) warped into the last image (light features) according to the induced homography H k l (on left) and the corresponding homography extracted from the super-homography H k (on right). Orange dashed lines are virtual lines that generate the SFPs in current (light '⋄') and last (dark '⋆') image.

Fig. 13 .

 13 Fig. 13. 2D trajectory of the vehicule motion along the Versailles sequence.

Fig. 14 .

 14 Fig. 14. Reconstruction of the road plane in the Versailles sequence.

Fig. 15 .

 15 Fig. 15.According to the hallway model, the modeling of urban canyons with three main planes allows a robust segmentation of the static scene. FPs extracted on areas which are not projections of these planes have to be viewed as obstacles.

Complete sequences are available on: http://caor.ensmp.fr/ ~simond

P k s ∝

Acknowledgment

This work was funded by the European project Cybercars. The authors thanks the IMARA (INRIA Rocquencourt) and LIVIC (INRETS-LCPC) projects which provide us the stereo-sequences we use.


V Jaumouillé

J-J Sinou

B Petitjean

Simulation of Payne effect of elastomeric isolators with a harmonic balance method

Keywords: Harmonic balance method, viscoelasticity, Payne effect

In the presented work, a non linear effect of rubber referred as Fletcher-Gent effect or Payne effect is investigated. It leads to a change in the rubber dynamic modulus with vibration amplitudes and, consequently, modifies resonance frequencies of mechanical systems including non linear elastomers. In this study a new methodology is developed to take into account Payne effect in a linear viscoelastic rubber material. Small vibration amplitudes around a no-preloaded state are predicted by considering frequency and amplitude dependencies of the material. This methodology has the advantage of using tabular experimental data from characterization tests which avoids the development of a complex model. In order to compute frequency responses, the non linear harmonic balance method is used and, for each iteration, new rubber properties are affected at each element according to its strain state. An equivalent strain measure is evaluated from the element strain energy density. This equivalent strain allows to associate dynamic properties of a material element subjected to multiaxial strain state with experimental dynamic properties of a material sample subjected to an uniaxial strain state. Practically, DMAP alter procedures are developed in order to evaluate energies in models defined with MSC.Nastran and the non linear solver is developed with Matlab. The method is applied on a satellite instrument isolator including four non linear rubber mounts. A non homogeneous spatial distribution of element equivalent strains is observed. Moreover, the maximum equivalent strain varies with frequency. These two observations validate the use of a specific methodology to deal with amplitude dependency of rubber.

Introduction

The use of rubber mounts is a common solution in industry in order to connect two structures and to achieve a combination of stiffness and vibration isolation. Kunz [START_REF] Kunz | Elastomer modelling for use in predicting helicopter lag damper behavior[END_REF] applied a rubber isolator technology on helicopters and Bruger et al. [START_REF] Bruger | FEM-Simulation of comfort relevant properties of elastorneric mountings: The Payne-effect under preloading[END_REF] optimized comfort properties of elastomeric properties for automotive industry. Simulations of rubber parts need to consider various complex effects and deal with properties varying with many parameters. Indeed, elastic and damping properties are strongly influenced by temperature and frequency (equivalence temperature-frequency), preload,... Many different experimental studies have been performed to study all these dependencies, more particularly on carbon-filled rubber which usually present a more complex behaviour [START_REF] Medalia | Effect of Carbon-Black on Dynamic Properties of Rubber Vulcanizates[END_REF][START_REF] Wang | Effect of polymer-filler and filler-filler interactions on dynamic properties of filled vulcanizates[END_REF]. One frequently studied dependency is the amplitude dependency. Mullins effect [START_REF] Mullins | Softening of rubber by deformation[END_REF] represents one of these effect where a softening effect is observed for the first excitation cycles on a rubber sample. A second effect, referred to Fletcher-Gent or Payne effect [START_REF] Payne | Low strain dynamic properties of filled rubbers[END_REF][START_REF] Fletcher | Non-linearity in the dynamic properties of vulcanisaed rubber compounds[END_REF], is the dependency of the dynamic properties on strain amplitude when material is submitted to small dynamic vibrations. Origin of Payne effect is not clearly established but supposed to be a continuous breaking and reforming Van der Waals forces between carbon-black aggregates [START_REF] Chazeau | Modulus recovery kinetics end other insights into the Payne effect for filled elastomers[END_REF]. Concerning the modelling of rubber, frequency dependent behaviour is today quite well modeled through the use of simple models involving combinations of springs and dashpots [START_REF] Ferry | ViscoelasticProperties of Polymers[END_REF][START_REF] Gil-Negrete | On the Modelling and Dynamic Stiffness Prediction of Rubber Isolators[END_REF]. For a full three dimensional case, hereditary integrals are used to formulate constitutive equations [START_REF] Lubliner | A Model of Rubber Viscoelasticity[END_REF]. Furthermore, fractional derivative models have also been developed to obtain a better frequency-dependent behaviour and to reduce the number of parameters to identify [START_REF] Enelund | Damping described by fading memory -analysis and application to fractional derivative models[END_REF]. Nevertheless, as regards the amplitude dependent behaviour modeling, more complex models have to be established for filled rubbers. The first one-dimensional model is the Kraus model [START_REF] Kraus | Mechanical losses in carbon-black filled rubbers[END_REF] who dealt with non linear viscous forces. Kraus considered contacts and shocks between aggregates. Later, Ulmer [START_REF] Ulmer | Strain dependence of dynamic mechanical properties of carbon black-filled rubber compounds[END_REF] improved this model to obtain a better approach of experimental observations. Other models tried to represent the irreversible slip processes occurring between filler particles by introducing friction elements [START_REF] Berg | A non-linear rubber spring model for rail vehicle dynamics analysis[END_REF][START_REF] Kaliske | Constitutive approach to rate-independent properties of filled elastomers[END_REF][START_REF] Rabkin | Material model and experimental testing of rubber components under cyclic deformation[END_REF]. These models are often limited to a one-dimensional case or are difficult to use in a frequency domain. Austrell et al. [START_REF] Austrell | A method to analyse the non-linear dynamic behaviour of carbon-black filled rubber components using standard fe codes[END_REF] made use of elastoplasticity elements to simplify implementation. The main drawback of these approaches is that a parameter identification is needed in conjunction with specific tests and that these models are limited to the time domain. For frequency analysis, these models are not very well adapted. However, some methodologies have been introduced by Gil-Negrete et al. [START_REF] Gil-Negrete | A simplified methodology to predict the dynamic stiffness of carbon-black filled rubber isolators using a finite element code[END_REF] who proposed to apply a first quasi-static step on a studied mount submitted to a given amplitude to assign adapted dynamic properties on each element according to its strain state in order to predict in a second time the dynamic stiffness of the mount at this same dynamic amplitude. Then, Rabkin [START_REF] Rabkin | Simulation of the Fletcher-Gent effect by using of the subroutine UPHI MSC/MARC (c)[END_REF] proposed to predict the dynamic stiffness of a rubber sample submitted to small vibrations around a preloaded state by modifying UPHI subroutine in a Marc analysis. This paper proposes to investigate this non linear Payne effect of rubber and develop a methodology to predict frequency responses of a complete system including non linear rubber mounts. This method does not include any parametric model and use tabular experimental data from a characterization test on a rubber sample. This approach uses a non linear harmonic balance method to deal with amplitude dependency. Only first harmonic is retained here due to practical observation [START_REF] Roland | Dynamic Mechanical-Behavior of Filled Rubber at Small Strains[END_REF]. However, more harmonics may be used in case of stronger non linearities. The technique is based on the calculation on an equivalent strain for each element and for each iteration in order to assign appropriate rubber dynamic properties. Indeed, each element of the rubber parts does not have the same strain state and consequently the same rubber property. This methodology is applied for models defined with MSC.Nastran and resolution is performed by the use of DMAP alter procedures and Matlab routines. Finally, it can be noted that the use of harmonic balance for problems in which the terms depend on the frequency and amplitude has been studied for the particular case of viscoelastic materials by Bilasse et al. : in [START_REF] Eh Boutyour | An approximated harmonic balance method for nonlinear vibration of viscoelastic structures[END_REF], the authors studied the nonlinear vibration of viscoelastic shell structures by coupling an approximated harmonic balance method with one mode Galerkin's procedure. In [START_REF] Bilasse | Linear and nonlinear vibrations analysis of viscoelastic sandwich beams[END_REF], they proposed numerical models with linear and nonlinear vibrations analyses of viscoelastic sandwich beams with various viscoelastic frequency dependent laws. The paper is organized as follows: method description is detailed in a first section. One application case is presented in a second section. In this application, a satellite instrument isolator including four rubber mounts is harmonically excited. Results are then shown to highlight the non linear behaviour of the system and the ability of the proposed technique to represent a non homogeneous spatial distribution of 2 element strains.

Methodology

This methodology proposes a new procedure to consider Payne effect for material modeling of viscoelastic materials in FE models with MSC.Nastran for small amplitude harmonic vibrations. Materials considered in this study are supposed to be isotropic, isothermal and nearly incompressible.

Viscoelasticity in MSC.Nastran for frequency analysis

For small amplitude strains, MSC.Nastran provides the ability to represent the behaviour of a viscoelastic material in the frame of the linear viscoelastic theory with a complex frequency-dependent shear modulus:

G(f ) = G ′ (f ) + iG"(f ) (1)
where G ′ and G" represent respectively the shear storage and loss moduli.

A viscoelastic behaviour leads to a complex dynamic stiffness matrix for all viscoelastic elements of the model. Stiffness matrix is computed on the basis of a representative reference modulus G ref which often corresponds to the long-term quasi-static shear modulus G 0 = lim ω→0 G(ω). G ref is independent of frequency. Thus complex stiffness matrix for viscoelastic elements is expressed as:

K = G ′ (f ) + iG"(f ) G ref K ref (2)
where K ref is the stiffness matrix with G ref as shear modulus. Additionally, a structural damping coefficient g may be affected on the whole model. A structural damping coefficient g ref must also be defined on the viscoelastic element group. To combine these two conditions, MSC.Nastran uses two tabular functions T R(f ) and T I(f ) to represent the complex moduli of all viscoelastic materials. The calculation of the complex stiffness matrix is then given by:

K = ((1 + g ref T R(f )) + i (g + g ref T I(f ))) K ref (3)
To assure equivalence with Equation (2), T R(f ) and T I(f ) are expressed as:

T R(f ) = 1 g ref G ′ (f ) G ref -1 (4)
T I(f ) = 1 g ref G"(f ) G ref -g (5)
MSC.Nastran does not provide the ability to define different dynamic properties for each particular viscoelastic element. When a viscoelastic material with amplitude dependent properties is subjected to a non homogeneous strain field, it may lead to severe errors for the dynamic stiffness matrix calculation.

Adaptation to Payne effect

The amplitude dependence of the complex shear modulus for small vibrations is referred to Payne effect and is noticeable for filled elastomers. This complex non linear effect is often attributed to irreversible slip processes between filler particles. Consequently, models of linear viscoelasticity are not suitable for frequency analysis. Based on an analogy with microplastic theory, some rheological models have been developed which introduce friction elements [START_REF] Berg | A non-linear rubber spring model for rail vehicle dynamics analysis[END_REF][START_REF] Kaliske | Constitutive approach to rate-independent properties of filled elastomers[END_REF]. These models are often formulated for uniaxial cases and are not adapted for multiaxial rubber models. Moreover a complex parameter identification has to be lead on material samples in terms of test conditions and parameter calculations. Up to now FE commercial codes do not provide the possibility to easily model Payne effect for frequency analysis.

An approximate procedure has been developed for dynamic analysis. This procedure has the advantage of using tabular experimental data from simple characterization tests (shear or compression tests) on material specimens and does not need any parameter identification for some simplified models. The non linear Harmonic Balance Method (HBM) is used to calculate frequency responses for a FE model integrating non linear viscoelastic elements (in terms of frequency and amplitude dependent shear modulus) subjected to an external dynamic load. Vibration amplitudes are not a priori known so equations have to be solved iteratively.

For each iteration of the solver, strain energy density is first computed element-wise and for the previous displacement field. Then, a local equivalent strain γ e eq is deduced so that the element strain energy density is equal to the homogeneous strain energy density in a characterization sample tested at the same dynamic amplitude γ e eq . This energy balance enables to associate dynamic properties of a material element subjected to multiaxial strain state with experimental dynamic properties of a material sample subjected to an uniaxial strain state. Afterwards, shear storage and loss moduli associated with this equivalent strain are deduced from experimental data and assembly of all elementary complex stiffness matrices is performed. A new displacement field is finally calculated from HBM equations and a non linear solver repeats this operation until convergence. In the following, HBM formulation is first presented in order to secondly detail the evaluation of non linear forces generated by all viscoelastic elements of the FE model.

HBM formulation

In order to compute frequency responses of a mechanical system integrating viscoelastic elements with amplitude dependent properties, HBM method is used. This method is able to compute periodic responses of non linear systems submitted to an external loading [START_REF] Jaumouillé | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems-Application to bolted structures[END_REF] or to self-excited vibrations [START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The constrained harmonic balance method, with application to disc brake squeal[END_REF]. HBM is also able to compute quasi-periodic responses [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] and able to conduct stability analysis [START_REF] Sinou | Effects of a crack on the stability of a non-linear rotor system[END_REF]. In the following, basic equations are recalled and a non linear solver is used for convergence. More details are given in Appendix 1 for the interested reader. Condensation techniques and continuation may also be applied in order to accelerate calculation. First we consider a discrete mechanical system with n degrees of freedoms (dofs) described with its n×n mass matrix M, stiffness matrix K and damping matrix D. The system includes non linear viscoelastic elements and v dofs are linked to these elements. Mass of all viscoelastic elements are integrated in the global mass matrix M. D and K include all linear parts of the model. Forces generated by the complex dynamic stiffness of all viscoelastic elements

K C v = K ′ v + iK ′′ v
are considered as an external non linear force F N L (X, Ẋ, Ω, t) which depends on degrees of freedom displacements X, velocities Ẋ, time t and angular frequency Ω. Ω is the angular frequency of an external periodic force F L (Ω, t) applied to the system. The governing equation of motion may be written as:

M Ẍ + D Ẋ + KX = F (X, Ẋ, Ω, t) = F L (Ω, t) + F N L (X, Ẋ, Ω, t) (6)
Then, we assume that the non-linear dynamical response of the system can be approximated by finite Fourier series:

X(t) = B 0 + m k=1 A k sin( k ν Ωt) + B k cos( k ν Ωt) X(t) = T(t)Z = I sin( Ω ν t)I cos( Ω ν t)I . . . sin( k ν Ωt)I cos( k ν Ωt)I . . . [B 0 A 1 B 1 . . . A k B k . . .] T (7)
where

I is the n × n identity matrix, Z = [B 0 A 1 B 1 . . . A k B k .
. .] T is the (2m + 1)n × 1 vector containing Fourier coefficients, m is the number of harmonics retained for the truncation, ν is an integer used to represent possible subharmonics, and

T(t) = [I sin( Ω ν t)I cos( Ω ν t)I . . . sin( k ν Ωt)I cos( k ν Ωt)I . . .]
is the n × (2m + 1)n matrix containing trigonometric functions. Finally, we assume that the vector force F N L (X, Ẋ, Ω, t) can be solved in finite Fourier series (see Appendix 1 for more details). Generally, the number of harmonic coefficients is selected on the basis of the number of significant harmonics expected in the non-linear dynamical response. Concerning Payne effect in viscoelastic materials, it has been described that, within certain strain amplitude limit, higher-order terms are not significant in the material response. The reader is referred to Roland [START_REF] Roland | Dynamic Mechanical-Behavior of Filled Rubber at Small Strains[END_REF] for more details. Consequently, a single harmonic has been retained in the HBM procedure leading to the following equation system (see Appendix 1 for more details):

K -Ω 2 M -ΩD ΩD K -Ω 2 M A 1 B 1 = b L (Ω) + b N L (Z v , Ω) (8)
where M, K and D are the mass, stiffness and damping matrices. b L (Ω) defines the Fourier coefficients from an external periodic force applied to the system. Ω is the angular frequency of the external periodic force. b N L (Z v , Ω) defines the Fourier coefficients from the non linear force F N L (X, Ẋ, Ω, t). It depends only on the Fourier coefficients Z v of the v degrees of freedom linked with the viscoelastic elements.

Non linear force evaluation

In order to calculate Fourier coefficients of non linear forces b N L , the non linear force F N L has to be detailed. Indeed, dofs are reorganized into v dofs which are linked to viscoelastic elements and nv dofs which are not using a boolean transition matrix P:

X = P X nv X v (9)
Thus, the complex dynamic stiffness matrix of all viscoelastic elements which is function of X v , may be expressed as:

K C v (X v , Ω) = K ′ v (X v , Ω) + iK ′′ v (X v , Ω) (10)
Thus, one may write:

F N L (X, Ẋ, Ω, t) = - 0 0 0 K ′ v (X v , Ω) X nv X v - 0 0 0 K ′′ v Ω (X v , Ω) Ẋnv Ẋv (11)
This stiffness matrix depends on the strain amplitude in the material and does not vary over one period. This remark allows to calculate non linear Fourier coefficients with the relation:

b N L = 1 T T 0 T T(t)F N L (X, Ẋ, Ω, t)dt b N L = -     0 0 0 0 0 K ′ v (Z v , Ω) 0 0 0 0 0 0 0 0 0 K ′ v (Z v , Ω)     Z -      0 0 0 0 0 K ′′ v Ω (Z v , Ω) 0 0 0 0 0 0 0 0 0 K ′′ v Ω (Z v , Ω)      ∇Z     S N Lnv S N Lv C N Lnv C N Lv     = -      0 0 0 0 0 K ′ v (Z v , Ω) 0 -K ′′ v Ω (Z v , Ω) 0 0 0 0 0 K ′′ v Ω (Z v , Ω) 0 K ′ v (Z v , Ω)          A 1nv A 1v B 1nv B 1v     (12)
Note that

Z v = [A 1v B 1v ]
T represents the Fourier coefficients of the v dofs in relation with all viscoelastic elements. Due to the amplitude-dependence of the stiffness matrix of all viscoelastic elements, a non linear solver is used. From a displacement field (i) Z obtained after iteration i, a new value of the viscoelastic dynamic stiffness matrix

(i) K v ( (i) Z v , f
) is performed in order to deduce a new value for the non linear forces

(i) b N L .
The strain field is not necessarily homogeneous in all the viscoelastic part leading to different material properties depending on the considered point. The new developed procedure is able to take into consideration this effect by affecting different material properties for each viscoelastic element. For an element e of the set V of all viscoelastic elements, an elementary matrix (i) K e v ( (i) Z e v , f ) is calculated from the displacement field (i) Z e v associated to this element. All these elementary matrices are then assembled to deduce the complete dynamic stiffness matrix (i) K v ( (i) Z v , f ):

(i) K v ( (i) Z v , f ) = e ∈ V (i) K e v ( (i) Z e v , f ) (13)
In the following, we detail the calculation of an elementary matrix (i) K e v ( (i) Z e v , f ) from an initial displacement field (i) Z e v .

Strain energy density

In the following, index i will be remove for simplification. All relations remain valid for each iteration.

For each viscoelastic element, an equivalent strain measure γ e eq is computed from an elementary displacement field Z e v . It allows to associate dynamic properties G ′e and G ′′e coming from an uniaxial test with the strain state in the element which combines different sollicitations. The strain equivalent strain is deduced at each iteration from the elementary displacement field Z e v and from an evaluation of the strain energy in an element calculated with a real shear modulus equal to G ref .

This strain energy may be expressed as:

E e (t) = 1 2 T X e v (t)K ref X e v (t) E e (t) = 1 2 T Z e v T T(t)T(t)N K ref Z e v E e (t) = 1 2 T Z e v sin 2 (Ωt)I sin(Ωt) cos(Ωt)I sin(Ωt) cos(Ωt)I cos 2 (Ωt)I N K ref Z e v E e (t) = 1 2 T Z e v 1-sin(2Ωt) 2 I sin(2Ωt) 2 I sin(Ωt) 2 I 1+cos(2Ωt) 2 I N K ref Z e v ( 14
)
E e (t) = 1 4 T Z e v N K ref Z e v + 1 4 T B e 1v K ref B e 1v -T A e 1v K ref A e 1v cos(2Ωt) + 2 T A e 1v K ref B e 1v sin(2Ωt)
Thus strain energy in the element is composed of a constant term

E e 0 = 1 4 T Z e v N K ref Z e v
which represents the mean value over one period and of an oscillatory term at the frequency 2Ω. The amplitude of the term is:

E e a = 1 4 ( T B e 1v K ref B e 1v -T A e 1v K ref A e 1v ) 2 + (2 T A e 1v K ref B e 1v ) 2 (15)
The peak value for each element over one period is consequently E e peak = E e 0 + E e a . The last step necessary for the evaluation of the equivalent strain is the computation of the strain energy density ρ e . ρ e is classically obtained by dividing E e peak by the element volume. Practically, MSC.Nastran is able to provide strain energy and strain energy density for each element of the model for a frequency analysis so strain energy density is obtained for all viscoelastic after introduction of a complex displacement field corresponding with Z v values at each iteration. An alter procedure in DMAP language is written to perform these operations.

Equivalent strain measure

Experimental data usually come from a simple shear or compression test and results in a table providing complex dynamic modulus depending on frequency and on a single dynamic strain amplitude γ exp . The strain energy density calculated at the previous step will allow to estimate an equivalent strain measure for each element γ e eq . This equivalent strain is comparable with the experimental strain amplitude and will allow to deduce appropriate storage and loss moduli to apply in each element. A energy balance is performed for each element by equalizing the obtained strain energy density with the homogeneous strain energy density in a characterization sample tested at a dynamic amplitude γ e eq . This amplitude is the equivalent strain amplitude which will be comparable with the experimental one γ exp . It has been shown in Gil-Negrete et al. [START_REF] Gil-Negrete | A simplified methodology to predict the dynamic stiffness of carbon-black filled rubber isolators using a finite element code[END_REF] that, for a material element subjected to its principal stresses, the deviatoric part of the total strain energy density may be expressed with the principal strains namely:

ρ d = G 3 (ǫ p1 -ǫ p2 ) 2 + (ǫ p2 -ǫ p3 ) 2 + (ǫ p3 -ǫ p1 ) 2 (16)
For our study, experimental data are associated with a compression test. Consequently, for a cylinder material element compressed with a strain amplitude of γ and with a shear modulus G ref , principal strains may be expressed as ǫ p1 = γ and ǫ p2,3 = -νǫ p1 = -γ 2 for a quasi incompressible material. In this case, distortion energy becomes:

ρ dT C = G ref 3 3 2 γ 2 + 0 + 3 2 γ 2 ρ dT C = 3G ref 2 γ 2 (17)
Energy density balance is then written for each element to deduce an equivalent strain measure:

ρ e dT C = ρ e 3G ref 2 γ e eq 2 = ρ e γ e eq = 2ρ e 3G ref ( 18
)
This procedure is an approximate procedure since the reference modulus is used to evaluate the equivalent strain. This choice has been made to keep a realistic value. Moreover, it has been shown with a virtual test on a material sample under compression that equivalent strain measure evaluated with the above procedure is equal to the homogeneous strain obtained in the material for this compression test. Then, a purely kinematic relation is not used as in Gil-Negrete et al. [START_REF] Gil-Negrete | A simplified methodology to predict the dynamic stiffness of carbon-black filled rubber isolators using a finite element code[END_REF] by convenience because MSC.Nastran provides the possibility to directly extract strain energy density of each element, avoiding to evaluate principal strains and averaging strain energy densities on all Gauss points. Finally, a dynamic calculation is performed to evaluate the equivalent strain instead of a quasi-static test in [START_REF] Gil-Negrete | A simplified methodology to predict the dynamic stiffness of carbon-black filled rubber isolators using a finite element code[END_REF] which leads to a better representativity on a larger frequency band.

Dynamic stiffness matrix assembly

From experimental tabular data, storage G ′e (γ e eq , f ) and loss G ′′e (γ e eq , f ) moduli are deduced using the equivalent strain measure and the frequency. Quadratic and cubic interpolations are employed to deduce approximate values for an equivalent strain and a frequency which do not exactly match experimental datas. This step is performed with Matlab. By analogy with Equation (2), each elementary matrix is calculated with the relation:

K e v (Z e v , f ) = K e v (γ e eq , f ) K e v (Z e v , f ) = G ′e (γ e eq , f ) + iG ′′e (γ e eq , f ) G ref K e ref ( 19
)
The final complex stiffness matrix of all viscoelastic elements K v (Z v , f ) is finally assembled and Fourier coefficients of the non linear force b N L is then computed according to Equation [START_REF] Enelund | Damping described by fading memory -analysis and application to fractional derivative models[END_REF]. For practical convenience, calculation of stiffness complex matrix and Fourier coefficients is performed with an alter procedure with MSC.Nastran which easily deals with connectivity of the model.

Practical procedure

The complete algorithm with the different continuation steps is presented in Figure 1. The algorithm is described in the particular case of a parameterized continuation with an curvilinear abscissa, a prediction with Lagrange polynomials and a Newton-Raphson correction method.

This procedure deals with the specific Payne effect in the "non linear term calculation". These detailed steps are then described in Figure 2. As seen previously, resolution has to be made incrementally. From an displacement field (i) Z v obtained at the iteration i, main steps to deduce the non linear force Fourier coefficients are synthesized here:

1. (i) Z v : displacement field of all viscoelastic elements obtained for iteration i.

2. Calculation of the strain energy density (i) ρ e (peak value) for each element. Shear modulus is equal to the reference modulus G ref .

3. Calculation of an equivalent strain measure (i) γ e eq = 2 (i) ρ e 3G ref with an energy balance. Evaluation of storage (i) G ′e ( (i) γ e eq , f ) and loss moduli (i) G ′′e ( (i) γ e eq , f ) with experimental data.

5. Calculation of each element complex stiffness matrix (i) K e v ( (i) Z e v , f ).

Assembly of the complex stiffness matrix of all viscoelastic element

(i) K v ( (i) Z v , f ).
7. Calculation of non linear force Fourier coefficients (i) A 1v and (i) B 1v .

It has also to be noted that this procedure is used for jacobian matrix calculation. Steps n • 1, 3, 4 are performed with Matlab whereas steps n • 2, 5, 6, 7 are performed with an alter DMAP procedure using with MSC.Nastran. 3 Application Case

Elastomeric isolator

An application is carried out on an elastomeric isolator from space industry. This system is designed to isolate a satellite instrument. Four rubber mounts are inserted between four metallic (aluminium) blocks clamped on the satellite and a circular part referred as "cross part" on which specific instrument is mounted. Figure 3 shows a detailed view of the considered structure. Satellite instrument is only represented here by a concentrated mass and does not appear on the figure. Weight of mounted instrument is approximatively 16kg. The global weight of the structure including isolator and instrument is nearly 20kg. Height of the structure is 180mm and cross diameter is 200mm. The four rubber mounts integrates two plane-parallel areas whom dimensions are 20 × 10 × 3.5 mm. Elastomeric parts are made of a rubber which presents a pronounced amplitude dependency. Material characterization has been performed for dynamic strain amplitudes of ±0.1%, ±0.25%, ±1%, ±2% and ±5%. Tested frequency bands are 1 -185Hz for ±5%, 1 -395Hz for ±2%, 1 -500Hz for ±0.1%, ±0.25% and ±1%. For simulations, only experimental data for a frequency band of 1 -185Hz will be considered in order to have material data for all dynamic amplitudes. In order to illustrate frequency and amplitude dependencies, experimental storage modulus is shown in Figure 4. Values have been removed for confidentiality. Isolator performance will be quantified with a force transmissibility. Indeed, for a given excitation force amplitude applied on the satellite instrument, reaction forces at the base of the four metallic blocks are observed. These blocks are clamped on the primary structure. For simulation, they are fixed. Force transmissibility is then deduced for different excitation configurations and maximum is observed for resonances. The main difficulty is the ability to predict resonance frequency and vibration amplitude due to the non linear stiffness of rubber mounts. This study focuses on the pumping mode so excitation is applied vertically on the satellite instrument. Consequently, rubber mounts are submitted to tension/compression.

Numerical model

A numerical model of the system is built with MSC.Nastran and is shown in Figure 3. This model is not initially submitted to a static preload. The satellite instrument is modelled as a single concentrated mass and an inertia matrix located at its center of gravity. This mass is rigidly connected to the upper surface of the cross part through rigid body elements which are not represented on Figure 3. Cross part and blocks are discretized with 10739 four-sided CTETRA and six-sided CHEXA solid elements. The geometry of elastomeric parts has been simplified in order to neglect elastomeric areas which are not significantly deformed. It has been checked that the change in the first six frequencies does not exceed 0.1%. Then a simplified mesh including only 192 CHEXA elements is created. A compromise has been established to keep a satisfying dynamic behaviour and a reasonable number of dofs. Change in the first six frequencies does not exceed 3% and the simplified mesh is sufficiently fine to avoid shear and volumetric mesh locking due to incompressibility. Consequently reduced integration elements or hybrid elements are not necessary. Moreover, condensation is applied on the cross part and on the four metallic blocks in order to reduce the number of dofs for simulation. First free mode of cross part and of each metallic block is higher than 1kHz and exceeds frequency band for simulation . This observation allows to carry out a Guyan condensation [START_REF] Rj Guyan | Reduction of Stiffness and Mass Matrices[END_REF] for these metallic parts on interfaces with rubber parts. After reduction, the model integrates 1662 dofs and 192 rubber elements. 6 dofs are used to represent the rigid body movement of the satellite instrument. Finally, an excitation force is applied on the single concentrated mass in order to obtain a significant non linear behaviour of the rubber mounts. A value of 145N is chosen in order to assure that the equivalent strain does not exceed 4% so that mechanical rubber properties remain in the range of experimental data. Indeed, for homogeneous rubber properties corresponding to an experimental dynamic amplitude of 5% (respectively 0.1%), a linear simulation with MSC.Nastran shows that the maximum equivalent strain obtained among all rubber elements (and for all the frequency band) is of 4% (respectively 2.41%). A value of 145N will consequently lead to a maximum equivalent strain measure included in this range of 2.41 -4% as desired.

Results

Simulations are performed on a frequency band of in order to capture the pumping mode.

In conjunction with the proposed non linear procedure, linear frequency analyses are carried out with MSC.Nastran with SOL108 algorithm. These calculations consider only frequency dependency on the viscoelastic material and affect homogeneous properties to rubber mounts corresponding to a fixed experimental dynamic amplitude.

Displacements

First results show that resonance of the pumping mode may vary on the frequency band of [75-100]Hz if homogeneous properties are considered. Indeed, Figure 5(a) describes the evolution of the displacement along z-axis of the concentrated mass for the six different dynamic amplitudes for linear simulations. Unsurprisingly, resonance frequency decreases for larger dynamic amplitude. The maximum value of 100Hz is obtained for an amplitude of 0.1% while the minimum of 75Hz is obtained for 5%. This decrease is associated with the rubber softening observed on experimental data in Figure 4. Maximum displacement also increases for larger dynamic amplitudes and varies between 67 and 111µm. Concerning the non linear procedure, resonance frequency is 83Hz which is an intermediary value showing that rubber is working in its non linear domain. The maximum displacement is 87µm and is coherent with linear calculations. Around resonance, this result is really close to the curve obtained with a homogeneous modulus chosen for an amplitude of 2%. However, it can also be noted that the curve shape is lightly different from linear calculations.

Transmissibility

A short analysis on the force transmissibility on the z-axis shows that the maximum value is around 3 for all simulations. The seven results vary between 2.96 and 3.37. This transmissibility is not constant and reaches a minimum value for a dynamic amplitude of nearly 2%. The non linear procedure is also very close to this value.

Equivalent strain

The proposed procedure is able to deal with non homogeneous properties of rubber and evaluates a different equivalent strain for each element allowing to affect different dynamic properties. This ability is illustrated on Figure 6(b) where a deformed view of a lower part of a rubber mount is shown. It corresponds to a compression phase. On each element, the equivalent strain is written so the spatial distribution is highlighted. Displayed colours do not correspond to a representation of the equivalent strain but to strain energy density of element surfaces. This view is a typical view from Patran. The maximum strain value of 3.21% is obtained for element 50099 and for a frequency of 83Hz, ie for resonance. As expected, this value is included in the range of 2.41 -4%. The rubber behaves in its non linear domain. The presented mesh is a simplified mesh and could be improved in order to increase the result precision. However, deformed shape shows that proposed mesh is sufficient to give a satisfactory behaviour of a nearly incompressible material and to avoid volumetric locking. Moreover, the equivalent strain is spatially distributed but also modified for each frequency. Figure 6(a) shows the evolution of the equivalent strain for the element 50099. The maximum value is obtained at resonance.

Conclusion

The main goal of this study is to propose a new approach in order to take into account Fletcher-Gent effect or Payne effect in a linear viscoelastic rubber material. It allows to perform dynamic simulations of systems integrating rubber parts to predict small vibration amplitude and resonances. No preload are considered here. The proposed methodology provides the ability to consider frequency and amplitude dependency by affecting a different material property to each different rubber material element. It allows to use tabular experimental data from characterization tests which avoids the development of a complex model. Due to the amplitude dependency, a non linear procedure called the harmonic balance method is used to calculate frequency responses. Consequently, an iterative procedure is applied at each frequency to determine vibration amplitudes which are not a priori known. In order to assign an adapted property to each element, an equivalent strain measure is evaluated from an energy balance. This energy balance enables to associate dynamic properties of a material element subjected to multiaxial strain state with experimental dynamic properties of a material sample subjected to an uniaxial strain state. Afterward, By replacing Eqn. [START_REF] Roland | Dynamic Mechanical-Behavior of Filled Rubber at Small Strains[END_REF] and Eqn. [START_REF] Jaumouillé | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems-Application to bolted structures[END_REF] into Eqn. [START_REF] Rabkin | Simulation of the Fletcher-Gent effect by using of the subroutine UPHI MSC/MARC (c)[END_REF], one obtains: MT(t)∇ 2 Z + DT(t)∇Z + KT(t)Z = T(t)b [START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The constrained harmonic balance method, with application to disc brake squeal[END_REF] Considering that for a n × n matrix W and a (2m + 1)n × 1 vector Y :

WT(t)Y = T(t)N W Y = T(t)diag(W, W, . . .) (2m + 1)n

×(2m + 1)n Y (26)
Equation ( 25) becomes:

T(t) N M ∇ 2 + N D ∇ + N K Z = T(t)b (27)
Time dependency may be suppressed and a frequency algebraic equation linking Fourier coefficient may be obtained using a Galerkin method which is a projection of the equation on trigonometric functions. Indeed these trigonometric functions define a scalar product:

< f, g >= 1 T T 0 f (t)g(t)dt (28)
Thus one may write:

1 T A may be expressed in a simpler manner:

A =              K K -Ω ν 2 M -Ω ν D Ω ν D K -Ω ν 2 M . . . K -k ν Ω 2 M -k ν ΩD k ν ΩD K -k ν Ω 2 M . . .              (31)
This system is equivalent of finding zeros of a function H : IR (2m+1)×n → IR (2m+1)×n :

H(Z) = A(Ω)Z -b(Z, Ω) = A(Ω)Z -b L (Ω) -b N L (Z v , Ω) (32)
We note that b is dependent on Z and Ω because b corresponds to the Fourier coefficients of F (X, Ẋ, Ω, t). b is a sum of Fourier coefficients from the linear excitation force b L and from non linear force b N L . b N L depends only on the Fourier coefficients Z v of the v dofs linked with the viscoelastic elements.

In the case where no analytical expression may be written between b N L and Z v , an evaluation of the approximate temporal terms X(t) and F N L (t) combining FFT procedures allows to numerically evaluate the Fourier coefficients b N L for each iteration. Concerning Payne effect in viscoelastic materials, it has been described that, within certain strain amplitude limit, higher-order terms are not significant in the material response [START_REF] Roland | Dynamic Mechanical-Behavior of Filled Rubber at Small Strains[END_REF]. Consequently, a single harmonic has been retained in the HBM procedure and ν = 1 leading to a 2n × 2n equation system:

K -Ω 2 M -ΩD ΩD K -Ω 2 M A 1 B 1 = b L (Ω) + b N L (Z v , Ω) (33)

Figure 1 :

 1 Figure 1: Complete HBM algorithm with non term procedure.

Figure 2 :

 2 Figure 2: Non linear force calculation for each iteration i

Figure 3 :

 3 Figure 3: Elastomeric isolator: (a) Global view, (b) Zoom on one rubber part. Rubber is coloured in red.

Figure 4 :

 4 Figure 4: Storage modulus: (a) Amplitude dependency; -(1Hz); --(45Hz); -.-(115Hz); ... (185Hz) (b) Frequency dependency; * (0.1%Hz); ⋄ (0.25%); o (0.5%); × (1%); ∇ (2%); + (5%).

Figure 5 :

 5 Figure 5: (a) Displacement of concentrated mass (b) Force Transmissibility along z-axis Rz F : -homogeneous properties of rubber without amplitude dependency);non linear procedure with amplitude dependency.

Figure 6 :

 6 Figure 6: Equivalent Deformation: (a) Equivalent strain for element id 50099 (b) Deformed rubber mount and equivalent strain distribution at 83Hz (figures correspond to equivalent strains and colours to strain energy densities).

  product on Eqn.[START_REF] Sinou | Effects of a crack on the stability of a non-linear rotor system[END_REF] leads to a (2m + 1) * n equation system:AZ = b with A = N M ∇ 2 + N D ∇ + N K (30)

shear storage and loss moduli associated with this equivalent strain are deduced from experimental data. Then, the non linear stiffness of rubber elements may be calculated at each iteration. The non linear solver repeats this operation until convergence. Practically, the methodology is developed to deal with models defined with MSC.Nastran. DMAP alter procedures are developed in order to evaluate energies in the model and non linear solver is developed with Matlab. The method is then applied on a satellite instrument isolator including four non linear rubber mounts. Results show a non homogeneous distribution of rubber properties since the equivalent strain differs from one element to another. Moreover, it has been shown that the maximum equivalent strain varies between 1.5% and 3.21% on a frequency band of [50-110]Hz. The Payne effect also strongly influences resonance frequency. Indeed, the pumping mode is located at nearly 83Hz while a procedure considering homogeneous properties may yield a resonance in the band [75 -110]Hz. Finally, the methodology can deal with any configurations of excitation (combination of torsion, shear...). In the case of a more pronounced non linearity, more than one harmonic may be chosen for simulations.

Appendix 1: HBM formulation

In this appendix, the HBM formulation is explained in details with application to the proposed study. We consider the previous governing equation of motion

where M, K and D are the mass, stiffness and damping matrices. Ω is the angular frequency of an external periodic force F L (Ω, t) applied to the system. F N L (X, Ẋ, Ω, t) concerns the nonlinear terms of the mechanical system. As previously explianed, the non-linear periodic response X(t) of the system can be approximated by finite Fourier series :

where

. .] T is the (2m + 1)n × 1 vector containing Fourier coefficients, m is the number of harmonics retained for the truncation, ν is an integer used to represent possible subharmonics, and

is the n × (2m + 1)n matrix containing trigonometric functions. The same work is then accomplished for the global force F :

In order to compute velocities and accelerations, a frequential derivative operator is defined:

Thus one may write: Ẋ(t) = T(t)∇Z and Ẍ(t) = T(t)∇ 2 Z (24)