J K Trømborg

H A Sveinsson

J Scheibert

K Thøgersen

D S Amundsen

A Malthe- Sørenssen

Slow slip and the transition from fast to slow fronts in the rupture of frictional interfaces

Keywords: friction, slow front, slow slip, multiscale modelling 1

The failure of the population of micro-junctions forming the frictional interface between two solids is central to fields ranging from biomechanics to seismology. This failure is mediated by the propagation along the interface of various types of rupture fronts, covering a wide range of velocities. Among them are so-called slow fronts, which are recently discovered fronts much slower than the materials' sound speeds. Despite intense modelling activity, the mechanisms underlying slow fronts remain elusive. Here, we introduce a multi-scale model capable of reproducing both the transition from fast to slow fronts in a single rupture event and the short-time slip dynamics observed in recent experiments. We identify slow slip immediately following the arrest of a fast front as a phenomenon sufficient for the front to propagate further at a much slower pace. Whether slow fronts are actually observed is controlled both by the interfacial stresses and by the width of the local distribution of forces among micro-junctions. Our results show that slow fronts are qualitatively different from faster fronts. Since the transition from fast to slow fronts is potentially as generic as slow slip, we anticipate that it might occur in the wide range of systems in which slow slip has been reported, including seismic faults.

Significance statement

Relative motion between solids in frictional contact is known to start progressively, with a slipping zone growing along the interface. The propagation of the front separating the stuck and slipping zones is usually very fast, but it can also slow down considerably. This transition

is not yet understood. Using a multi-scale model, we demonstrate how the transition to slow fronts is the direct consequence of slow slip motion at the interface, and that slow fronts are qualitatively different from faster fronts. Since the transition to slow fronts is allowed by slow slip, we expect that it can occur in the wide range of systems in which slow slip has been reported, from engineering to seismology.

Introduction

The rupture of frictional interfaces is a central mechanism in many processes including snow slab avalanches, human object grasping and earthquake dynamics [START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF]. Rupture occurs through the propagation of a crack-like micro-slip front -the rupture front -across the interface. This front represents the moving boundary between a stick region and a slipping region that coexist within the contact. In so-called partial-slip situations, fronts propagate quasi-statically at a pace controlled by the external loading, as studied in mechanical engineering for decades [START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Prevost | Probing the micromechanics of a multicontact interface at the onset of frictional sliding[END_REF]. Recently, fast cameras enabled the observation of much faster fronts, which are classified into three types: supershear fronts faster than the material's shear wave speed c s , sub-Rayleigh fronts propagating at velocities close to c s and slow fronts much slower than c s [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF][START_REF] Xia | Laboratory earthquakes: the sub-Rayleighto-supershear rupture transition[END_REF][START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF][START_REF] Nielsen | Experimental observation of stick-slip instability fronts[END_REF][START_REF] Audry | Slip dynamics at a patterned rubber/glass interface during stick-slip motions[END_REF]. Whereas the two first types have been predicted theoretically, the physical mechanisms underlying slow fronts are still debated.

A better understanding of slow fronts appears as a significant step towards an improved assessment of how frictional motion begins. It is also expected to shed light on the important topic of slow earthquakes, which have been increasingly reported in the last decade [START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF]. In this context, an intense theoretical and numerical activity arose to investigate the origins and properties of rupture fronts. Two different approaches have been explored.

On the one hand, two-or three-dimensional elastodynamic models have been used to relate the macroscopic loading conditions to the stress field along the contact interface [START_REF] Ben-Zion | Dynamic ruptures in recent models of earthquake faults[END_REF][START_REF] Bartolomeo | Dynamic rupture at a frictional interface between dissimilar materials with asperities[END_REF][START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF][START_REF] Kammer | On the propagation of slip fronts at frictional interfaces[END_REF][START_REF] Otsuki | Systematic breakdown of Amontons' law of friction for an elastic object locally obeying Amontons[END_REF]. These local stresses were indeed shown experimentally to play a role in the selection of the front type [START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF]. However, the models were based on simple friction laws, e.g. Coulomb friction [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF] or velocity-weakening friction [START_REF] Kammer | On the propagation of slip fronts at frictional interfaces[END_REF][START_REF] Otsuki | Systematic breakdown of Amontons' law of friction for an elastic object locally obeying Amontons[END_REF] and did not exhibit slow fronts. On the other hand, fronts much slower than the speed of sound were produced by one-dimensional friction models using improved local friction laws involving an intrinsic time-scale: either an aging scale (i) related to the long-time strengthening of an interface at rest [START_REF] Sinai | Slow rupture of frictional interfaces[END_REF][START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF] as is classically considered in rate-and-state theories [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF][START_REF] Marone | Laboratory-derived friction laws and their application to seismic faulting[END_REF], or a dynamic scale (ii) related to the short-time collective dynamics of formation and rupture of a statistical number of micro-junctions during the rupture of the interface [START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Capozza | Static friction and the dynamics of interfacial rupture[END_REF]. These two time-scales (i) and (ii) were shown, in the same experiment [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF], to control the recovery of contact area after slip arrest and the transition between a fast and a slow regime for the slip motion triggered by the front passage, respectively. However, because one-dimensional models are unable to reproduce realistic stress distributions at the interface, they do not allow for quantitative comparison with experiments. Here we combine both 2D elastodynamics and time-dependent friction into a multi-scale model for rupture fronts. We demonstrate that the model simultaneously reproduces two separate, unexplained experimental observations: the transition between fast and slow front propagation during a single rupture event; and the transition from fast to slow slip motion at the interface shortly after rupture. In addition, through a more general study of the model, we show that slow slip occurring immediately after a fast front's arrest is a sufficient phenomenon for the front to propagate further as a slow front, and we identify parameters involved in front type selection.

Model description

The frictional stability of a system made of two solids in contact depends on the level of normal and shear stresses at the contact interface. These interfacial stresses result from the external forces applied at the boundaries of the solids, transmitted through the bulk. Slip motion will in general be triggered when the local interfacial shear stress reaches a threshold, the level of which crucially depends on the interface behavior law at the microscale.

The net contact between two solids generically consists of a large number of stress bearing micro-junctions whose nature depends on the type of interface. For rough solids, each microjunction corresponds to a micro-contact between antagonist asperities, whereas for smoother surfaces the junctions can be solidified patches of an adsorbate layer [START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF]. The three physical aspects of the junction behavior that we consider to be essential are as follows. 1) A microjunction in its pinned state can bear a shear force f T , provided it remains smaller than a threshold f thres . When f thres is reached, a local fracture-like event occurs, and the junction enters a slipping state. 2) In the slipping state, the micro-junction can let the interface slip, either through the micro-slipping of micro-asperities in contact or through the fluidization of an adsorbate layer. During slip, the micro-junction sustains some residual force f T = f slip , with f slip smaller than f thres . 3) Slipping micro-junctions have a certain probability to disappear or relax.

For example, a micro-contact disappears when an asperity moves away from its antagonist asperity by a typical distance equal to the mean size of micro-contacts, as classically considered for slow frictional sliding, e.g. in rate-and-state friction laws. However, another picture may arise as a consequence of the sudden release of energy when pinned junctions break. This energy will transiently heat the region around the micro-junction [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF]. The rise in temperature will significantly increase the rate of a thermally activated relaxation of the slipping micro-junction during the time necessary for the interface to cool down [START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF]. The frictional consequences of such temperature rises have recently received renewed attention (see e.g. ( 22)), but remain poorly understood. In an attempt to include such thermal processes in our model, we recognize that they will lead to time-rather than distance-controlled relaxations, so that the shear force drop will be distributed in time. In order for the interface to continue bearing the normal forces applied to it, the micro-junctions that relax are replaced by new, pinned junctions bearing a small tangential force f new .

The physical aspects described above have been modeled in a simple way using the following assumptions (Methods). We consider the rough frictional interface between a rigid track and a thin linear elastic slider of length L and height H (Fig. 1A). The bulk elastodynamics of the slider are solved using a square lattice of blocks connected by internal springs [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF][START_REF] Yim | Numerical simulation and visualization of elastic waves using mass-spring lattice model[END_REF] (Fig. 1B). The multi-contact nature of the interface is modeled through an array of N s tangential springs representing individual micro-junctions, attached in parallel to each interfacial block (Fig. 1C) [START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Capozza | Static friction and the dynamics of interfacial rupture[END_REF][START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF]. The individual spring behavior is as follows (Fig. 1F, [START_REF] Thøgersen | History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework[END_REF]). A spring pinned to the track stretches linearly elastically as the block moves, acting with a tangential force f T on the block. When the force reaches the static friction threshold f thres (we neglect aging, so that f thres is time independent), the micro-junction ruptures and the spring becomes a slipping spring acting with a dynamic friction force * f T =f slip . After a random time t R drawn from a distribution Τ(t R ), the slipping spring relaxes. It is replaced immediately by a pinned, unloaded spring (f new =0) and a new cycle starts. Here we use Τ(t R ) as a simplified way of modeling the distribution of times after which micro-junctions relax. Due to the variety and the complexity of the underlying thermal processes, we did not try to derive Τ(t R ) for a specific situation. Rather, we chose to model Τ(t R ) in the simplest way, as a Gaussian with average time <t R > and width δt R . The shape of Τ(t R ) is not crucial: we obtain qualitatively similar results with an exponential distribution. The width of Τ(t R ) is the only source of randomness in our model and causes the interface springs of a block to evolve differently from each other.

The model reproduces rupture front observations

In this section we use the model to reproduce unexplained experimental observations made in polymethylmethacrylate by Fineberg and coworkers. We use the loading conditions and material constants of (4), as done in [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF]. We then select the parameters of the friction law in order to reproduce both the complex space-time evolution of the front velocity reported in [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF] and the time-dependent interfacial slip history reported in [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF].

* f thres and f slip are taken proportional to the normal force p on the corresponding block (Methods). This assumption ascribes the pressure dependence of the forces on a block to the individual forces bore by a constant number N s of springs, rather than to a pressure dependent number of springs per block.

The slider is first submitted to normal load F N only. Then, a slowly increasing tangential load F T is applied on its left side, at height h above the interface (Supporting Methods).

Macroscopically (Fig. 1D), the slider is first loaded elastically, yielding high shear stresses in the vicinity of the loading point. These growing stresses eventually trigger the slip of a first block, the motion of which increases the force on its neighbors (Fig. 1E) until they also start to slip, and so on. This mechanism results in the propagation of a rupture front across the interface. The first such events arrest before reaching the contact's leading edge and correspond to so-called precursors to sliding [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF][START_REF] Rubinstein | Dynamics of precursors to frictional sliding[END_REF][START_REF] Scheibert | Role of friction-induced torque in stick-slip motion[END_REF][START_REF] Maegawa | Precursors of global slip in a longitudinal line contact under non-uniform normal loading[END_REF][START_REF] Amundsen | 1D model of precursors to frictional stick-slip motion allowing for robust comparison with experiments[END_REF]. In the following, we will only consider interface-spanning events (Fig. 2), that are responsible for the large force drops in Fig. 1D. We measure the front propagation speed v c as a function of position x along the interface (Fig. 2C) by defining rupture at the block scale, when the fraction of pinned springs drops below a threshold value (Figs. 2A,B).

Figure 2A shows an event starting as a fast front (v c ~cs /3), then turning into a slow front (v c ~cs /100) at x~0.4L, before turning back into a fast front at x~0.6L. This space-time development is in excellent agreement with the experiments from which the model parameters were taken (4). In particular, the locations of the transitions between fast and slow fronts, the duration of the slow front (~1ms) and the velocity ratio between fast and slow fronts are all matching the experimental observations. Other events are fast across the whole interface (Fig. 2B). We emphasize that such spatially heterogeneous dynamics arise spontaneously in our model. This success is presumably due to the two-dimensionality of the model, which allows it to reproduce interfacial stress heterogeneities arising from macroscopic sample geometry and external loading (11) (Fig. S1).

Figure 3C shows the slip experienced by an interfacial block during an event, the timeevolution of which exhibits two distinct regimes. A fast slip regime is followed by a roughly linear slow slip regime, in good agreement with the slip history reported in [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF]. Notably, both the ratios of slip distances and of slip velocities experienced in either regime are matching the experimental observations. The simultaneous agreement found with two independent measurements -reproducing the fronts' spatio-temporal dynamics at the slider scale and the slip dynamics at the block scaleindicates that our model can provide insight into the set of experiments reported by Fineberg and coworkers. Because the basic ingredients of the model represent general features of frictional interfaces, we expect the model to be applicable to a larger class of systems.

Relationship between slow slip and the transition to slow fronts in the model

From this section on we change the scope from reproducing particular experiments to performing a systematic study of the model. We begin by varying select model parameters to unravel the mechanisms underlying slow fronts.

Figure 3C (reference curve) shows the slip history of a block during an event with a slow front (Fig. 3A); it has both fast and slow slip regimes. Fast slip initiates with the passage of the rupture front and is independent of <t R > (the three curves in Fig. 3C overlap during fast slip).

It is followed by a short period of arrest. Then, slow slip originates from the following mechanism: For each spring leaving the slipping state after t R the friction force on the block is reduced. As a result the block moves a small distance, driven by the forces from the neighboring blocks. Slow slip results from the succession of such incremental block movements. A scaling analysis of this mechanism (Supporting Equations) yields a prediction for the slow slip velocity v slowslip , which quantitatively captures the data (inset of Fig. 3C).

The fast-slow-fast front in (A) becomes an arrested front (B) when we restart the simulation in the same state, but with slow slip turned off. We do this by setting f new = f slip so that springs leaving the slipping state no longer relax, but return to the pinned state bearing the same force as during slip. The front in (B) arrests where the fast-to-slow transition used to take place. These results support the following scenario: Slow fronts are fronts that would arrest in the absence of slow slip, but are pushed further by the increasing loading due to slow slip behind the front tip.

To show that slow slip occurring immediately after front arrest is a sufficient condition for the transitions to slow fronts to be allowed, we consider the following simplified model. We leave microscale dynamics out by replacing the ensemble of springs connecting each block to the track with a single spring: The spring breaks when reaching the force N s f thres , it slides with a force N s f slip and it returns to the pinned state still bearing the force N s f slip when the block velocity vanishes. Slow slip is introduced directly as a small velocity applied to the spring attachment point (Supporting Methods). The transition from fast to slow fronts is indeed observed in this simplified model (Fig. S2), and disappears if slow slip is turned off.

To characterize the relationship between slow slip velocity and slow front velocity, we consider both the complete and simplified models and use interfaces prepared as follows (Supporting Methods). We choose the normal force p(x) on each block from the substrate to be uniform in space and constant in time, and the shear force profile τ(x) to have initial values τ0(x) uniform in space except in a region used to initiate front propagation (Fig. S3). We vary different model parameters, keeping the front initiation region unchanged, and systematically measure the slow front velocity v c,slow as a function of v slowslip (Fig. 3D,inset). Suitable rescaling (Supporting Equations) of v slowslip allows all points to collapse on a single straight line (Fig. 3D). This collapse shows that v c,slow is proportional to v slowslip , with the same conversion factor for both models. We conclude that even though v slowslip depends on the interfacial dynamics, the conversion factor does not; it mainly depends on the way stresses are transferred from the slowly slipping region to the rupture tip. Similar linear relationships between slip and front velocities were found in previous models [START_REF] Sinai | Slow rupture of frictional interfaces[END_REF][START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF][START_REF] Ohnaka | A cohesive zone model for dynamic shear faulting based on experimentally inferred constitutive relation and strong motion source parameters[END_REF].

Selection of the front type in the model

We now turn to the question of front type selection. Ben-David et al. [START_REF] Ben-David | The dynamics of the onset of frictional slip[END_REF] showed that slow (fast) fronts correlate locally with small (high) shear to normal stress ratio τ0(x)/p(x) just before the event. In equilibrium, τ is the sum of the spring forces f T for each block. Due to the randomness in the slip times t R the individual f T will be different. As in previous works [START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF][START_REF] Thøgersen | History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF] it is useful to define a distribution φ(f T ) of the forces in the springs attached to each block. φ is a dynamical quantity evolving with both the loading on and the motion of blocks.

Note that in general φ will be different for every block.

We find that the width σ of the distribution is as important as the stress ratio for front type selection. To demonstrate this we take event II as a reference and modify it by changing σ for all blocks at x>0.3L. To isolate the effect of σ from the effect of the stress ratio on front propagation we increase σ, but keep τ0(x) unchanged. We do this by making each block's φ(f T ) a uniform distribution while keeping its average unchanged (Fig. S4). We then restart the simulation and observe that the front has changed from a fast-slow-fast one (Fig. 4A) to a fast one across the whole interface (Fig. 4B). The reason for this behavior is that increasing σ from Fig. 4A to Fig. 4B weakens the interface, which enables the fast rupture to reach the leading edge without stopping and becoming a slow front. The weakening can be understood from Fig. 4C, which shows the effective static friction coefficient µ s eff =τ max /p of an interface block, with τ max the maximum friction force on the block before it starts to slip. Small (large)

σ correspond to large (small) τ max (Fig. S5), i.e. a locally stronger (weaker) interface, consistent with the fact that points corresponding to slow fronts cluster at low σ in Fig. 4C.

Figures 4A and4B indicate that, besides τ0/p, σ is involved in front type selection. Using interfaces prepared as in Fig. 3D, we systematically vary the values of τ0 and σ and observe which front type is selected (Fig. 4D). The effect of τ0 agrees with experiments ( 6), but is modulated by the effect of σ. Low τ0 and σ lead to fronts that arrest before spanning the whole interface. Large τ0 and σ yield fast fronts. Global events containing a slow part are found in a region of intermediate τ0 and σ. These results are found qualitatively robust not only against changes in sample geometry or energy stored in the nucleation region, but also against changes in the individual behavior of the interfacial springs (Fig. S6). We therefore expect the main features of Fig. 4D to be widely relevant.

Discussion

The present model differs from classical rate-and-state friction laws [START_REF] Baumberger | Solid friction from stick-slip down to pinning and aging[END_REF][START_REF] Marone | Laboratory-derived friction laws and their application to seismic faulting[END_REF], which are empirical laws based on experimental results obtained for small sliding velocities (≲100µm/s), for which self-heating of the interface is negligible. In rate-and-state models the relevant time-scale is that of aging, i.e. the slow (logarithmic) recovery of contact area at rest.

The relevant slip distance is the average micro-contact size, after which all junctions are renewed. Here, we focus on a different friction regime, with a µs-long fracture process, followed by large slip velocities in the fast slip regime (~100mm/s in (20)), both phenomena contributing to a significant temperature increase at the interface. Our model thus considers the relevant time-scale t R to be the one of thermally activated relaxations during the period needed for the interface to slowly slip and stop (we neglect aging). The corresponding slip distance is the one allowed by junction relaxations during the same period.

In the context of rate-and-state friction, it was recently suggested that slow slip velocity corresponds to the sliding speed at which the interface's steady state friction laws reaches a minimum [START_REF] Sinai | Slow rupture of frictional interfaces[END_REF][START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF][START_REF] Kaproth | Slow earthquakes, preseismic velocity changes, and the origin of slow frictional stick-slip[END_REF]. In our model, however, slow slip has a completely different origin, related to the intrinsic relaxation dynamics of the interface after a fast slip period (Fig. 3C, Supporting Equations). As a matter of fact, the steady-state friction law that emerges from the micro-junction dynamics at the block scale is purely monotonic and decreasing (Fig. 1G, [START_REF] Thøgersen | History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework[END_REF]).

To what extent the slow fronts related to both types of slow slip share similar properties remains to be investigated.

Up to now, the term slow front has been used in the literature to name any front propagating at least one order of magnitude slower than the material's wave speeds. Our results allow for a tentative definition of the slow fronts first observed in [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF]. They are fronts that would arrest in the absence of slow slip, but are pushed further by the increasing loading due to slow slip behind the front tip. Importantly, the motion results from the intrinsic relaxation dynamics of the interface and not from an external driving. A direct consequence is shown in Fig. 5A, where the increase in tangential load is stopped during the propagation of the slow front of Fig. 2. We observe that the slow front is not affected and propagates exactly like in the original event, which indicates that the front is truly dynamic. This contrasts with the behavior of the event shown in Fig. 5B, where the loading is applied to the bottom-left-most block.

Extremely slow fronts are then observed, similar to the ones reported in e.g. [START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF][START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Capozza | Static friction and the dynamics of interfacial rupture[END_REF][START_REF] Bar-Sinai | Instabilities at frictional interfaces: Creep patches, nucleation, and rupture fronts[END_REF], which propagate over the time-scale of the interval between two dynamic events. These fronts stop when the driving stops, indicating that they are of a different, quasi-static nature.

This distinction sheds light on the difference between our model and recent models for the onset of sliding of extended interfaces. Urbakh and coworkers [START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Capozza | Static friction and the dynamics of interfacial rupture[END_REF] have used 1D springblock models with an interfacial behavior also based on micro-junctions having two possible states. However they used f slip =0, which prevented slow slip to occur. As a consequence, they observed fast and quasi-static fronts, but did not observe dynamic slow fronts like those reported here. Bouchbinder and coworkers [START_REF] Sinai | Slow rupture of frictional interfaces[END_REF][START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF][START_REF] Bar-Sinai | Instabilities at frictional interfaces: Creep patches, nucleation, and rupture fronts[END_REF] have developed an improved rate-andstate law and used it in 1D. They consider the classical aging time scale whereas we consider a dynamic healing time scale relevant during the short-lived rupture-induced temperature rise of the interface. They observed the transition from quasi-static to fast dynamic rupture, but no fast to slow dynamic front transition. Other works considered 2D models, but with velocity-weakening local friction laws [START_REF] Kammer | On the propagation of slip fronts at frictional interfaces[END_REF][START_REF] Otsuki | Systematic breakdown of Amontons' law of friction for an elastic object locally obeying Amontons[END_REF]. While they observed fast and quasi-static fronts, they did not report any transition from fast to slow dynamic fronts.

Our results make a direct bridge between the separate observations of a transition from fast to slow fronts and of slow slip in recent experiments [START_REF] Rubinstein | Detachment fronts and the onset of dynamic friction[END_REF][START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF]. As slow slip was observed in many systems from geoscience [START_REF] Peng | An integrated perspective of the continuum between earthquakes and slow-slip phenomena[END_REF][START_REF] Hirose | Repeating short-and long-term slow slip events with deep tremor activity around the Bungo channel region, southwest Japan[END_REF][START_REF] Ohnaka | Characteristic features of local breakdown near a crack-tip in the transition zone from nucleation to unstable rupture during stick-slip shear failure[END_REF] to materials science [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF][START_REF] Heslot | Creep, stick-slip, and dry-friction dynamics: Experiments and a heuristic model[END_REF], we expect the transition to slow fronts to be possible in these systems too. The physical process underlying slow slip, e.g.

thermal softening in glassy polymers [START_REF] Ben-David | Slip-stick and the evolution of frictional strength[END_REF], thermal creep in paper [START_REF] Heslot | Creep, stick-slip, and dry-friction dynamics: Experiments and a heuristic model[END_REF] and dilatant strengthening in subduction zones [START_REF] Segall | Dilatant strengthening as a mechanism for slow slip events[END_REF] differs from system to system [START_REF] Toro | Fault lubrication during earthquakes[END_REF]. The present model does not aim at modeling one particular process. Rather, through the time distribution T(t R ),

we introduce the minimal generic ingredient necessary to yield a slow interfacial slip motion.

Using the simplified model, we demonstrated that slow fronts can in principle exist even if only one micro-junction is kept per block, i.e. the force distribution φ has width σ=0, provided slow slip is introduced in another way. However, we believe that virtually all frictional interfaces are of the multi-junction type at the relevant scale, e.g. micro-contact or molecular scale. Due to random physical properties or stochastic pinning/relaxation events combined with the previous sliding dynamics there will always be some disorder in the forces bore by the junctions, as recognized in various friction models (e.g. [START_REF] Braun | Dynamics of transition from static to kinetic friction[END_REF][START_REF] Persson | Theory of friction: Stress domains, relaxation, and creep[END_REF][START_REF] Thøgersen | History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework[END_REF][START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]). The disorder was shown to control e.g. the steady sliding friction force and the transition between smooth sliding and stick-slip regimes [START_REF] Braun | Master equation approach to friction at the mesoscale[END_REF]. Here we showed that the force distribution is also a key parameter in front type selection. It appears as a state parameter that, by controlling the effective friction threshold µ s eff at the block scale [START_REF] Thøgersen | History-dependent friction and slow slip from time-dependent microscopic junction laws studied in a statistical framework[END_REF], is able to affect interfacial rupture in a way analogous to but different from the usual age state used in rate-and-state friction laws.

We believe that all our results are relevant to a large class of systems with a random population of stress-bearing entities having both a threshold like rupture behaviour and a time-dependent healing process, like in nanoscience [START_REF] Li | Frictional ageing from interfacial bonding and the origins of rate and state friction[END_REF], polymer science [START_REF] Cordier | Self-healing and thermoreversible rubber from supramolecular assembly[END_REF] and seismology [START_REF] Mclaskey | Fault healing promotes high-frequency earthquakes in laboratory experiments and on natural faults[END_REF].

Methods

Parameters are in Supporting Table S1. Spring j of block 𝑖 has strength 𝑓 !!!"#! 𝜇 ! 𝑝 ! /𝑁 ! and stiffness 𝑘 !" ∝ 𝑝 ! /𝑁 ! ,. The force on block 𝑖 from spring 𝑗 is 𝑓 !"# = 𝑘 !" (𝑥 !" -𝑥 ! ), where 𝑥 !" is the attachment point of the spring to the track and 𝑥 ! is the block's position. While slipping, the spring trails the block with 𝑥 !" adjusted in every time step to ensure 𝑓 !"# = 𝑓 !"#$ = 𝜇 ! 𝑝 ! /𝑁 ! . Note that if the block motion reverses, 𝑥 !" stops changing and remains fixed while

𝑓 !"# ≤ 𝑓 !"#$ . The probability distribution 𝑇 𝑡 ! is based on a Gaussian 𝑇 ! 𝑡 ! = ! !!!! ! exp (- (! ! !! ! ) ! !!" ! ! )
, modified so that negative lifetimes are forbidden:

𝑇(𝑡 ! ) = 𝑇 ! 𝑡 ! + 𝑇 ! -𝑡 ! , 𝑡 ! ∈ [0, ∞].
With the parameter values we use, < 𝑡 ! >≅ 𝑡 ! .

In the simulations behind Fig. 3C and 3D, the ratio 𝛿𝑡 ! /𝑡 ! is maintained. Details on initialization and non-frictional boundary conditions are in Supporting Methods. The motion is found using a Velocity Verlet scheme with Δ𝑡 = 2 • 10 !! s.

μ eff s = ððτ max =pÞ -ðf slip =f N ÞÞ= ððf thres =f N Þ -ðf slip =f N ÞÞ vs. σ = σ=ðf thres -f slip Þ
for all blocks in events I (green) and II (red), and for prepared events (black) with homogeneous τ 0 = ððτ 0 =pÞðf slip =f N ÞÞ=ððf thres =f N Þðf slip =f N ÞÞ. •, slow front; ×, fast front; ○, loading zone (x < 2 cm). Solid (dashed) line is for uniform (bell-shaped) ϕ(f T ) (Fig. S4). (D) Observed front type for prepared interfaces vs. both τ 0 and σ. Dots represents data points.

Supporting Methods

This section supplements the model description found in the main text and its Methods section with detailed information on how we initialize the system and apply the boundary conditions. We also provide some additional information on the simplified model.

The slider is initialized with full normal load F N and no tangential load F T by gradually applying F N without allowing springs to break, a technicality required because the normal forces on the springs, f N ij , start at zero and therefore springs, if allowed to, would break under any stretching. We distribute the load F N uniformly on the top blocks; apart from this we use the same non-frictional boundary conditions as in [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF]. The unique equilibrium is found through damped relaxation of typical duration 10 ms. After relaxation, we check that no spring is stretched beyond its strength and introduce the driving spring starting from zero applied driving force F T . Then F T , which acts on the block on the left side of the slider situated at height h above the interface, through the driving spring, increases as the driving point moves to the right with speed V .

In the simplified model used for Fig. 3D, we disregard the microscopic state by using a single friction spring per block. Taking parameters from the microscopic reference model described in the main text, each block's spring now has a strength τ thres = µ s p i = N s f thres . The stiffness k i of a block's friction spring equals the combined stiffness of the springs per block in the reference model. The force on block i from its friction spring is

f T i = k i (x is -x i )
, where x is is the attachment point of the spring to the track. Upon breaking, the spring becomes a slipping spring and its behaviour starts to differ from that of the springs in the microscopic model. We impose a slow slip by letting x is move with a velocity v slow slip spring for a time t slow slip spring = tR . This process competes with the dynamic friction law where the spring trails the block with x is adjusted in every time step to ensure |f T i | ≤ µ d p i , so that the spring attachment point moves with the highest of v slow slip spring and v xi , the speed of the block in the x-direction. When the block motion reverses (v xi changes sign) the spring returns to the pinned state, but x is continues to move at v slow slip spring until t slow slip spring after the spring broke.

The systematic studies leading to Fig. 3C (inset, green), 3D and 4D were done with different normal forces and different initialization from the other simulations. The normal force boundary conditions on the top and bottom were exchanged: this simplifies the analysis by setting a constant normal force p i = F N /N x on all blocks i at the interface. To maintain stability against global rotation, the top blocks interacted with an elastic ceiling with the same properties as the elastic foundation used in Ref. [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF] and the other simulations presented here.

To obtain an initial state with a prescribed interfacial shear stress profile we turned the interface springs off during the initialisation. In their place we added to each bottom block the force corresponding to the shear stress to be prescribed. We also introduced the driving spring, but let V = 0. During relaxation, the sample moved along the x-axis until the force in the driving spring balanced the net force from the interfacial shear stress. To get rid of oscillations more efficiently we added damping forces -α( v i ) on the blocks' motion. After relaxation, the extra forces and the extra damping were turned off and the interfacial springs were introduced, with their attachment points x ij chosen such that the net force on each block was unchanged and the desired distribution of spring forces, φ(f T ), appeared. We then waited a few timesteps to ensure that the transition from pre-to post-relaxation involved no force discontinuities. Next, instead of driving the system with V = 0 until rupture was triggered, we started fronts by depinning simultaneously all springs for all blocks to the left of x trigger . The shear stress in the triggering region has a strong influence on the rupture fronts; in order to compare results between simulations we used a constant value τtrigger . Supporting Figure S1: Spatial distributions of normal (p) and shear (τ ) stresses from the simulation behind events I-III. (A) The initial state with full normal load F N and no tangential load F T (see Supporting Methods). The normal stress p(x) is symmetric, with edge effects related to the flat punch geometry of the contact. The shear stress τ (x) is antisymmetric, due to Poisson expansion being restricted at the interface by friction. (B) The state just before event I. The application of F T has modified both the shear stress profile and the normal stress profile (due to the friction-induced torque arising when F T is applied at a finite height h above the interface). In order to isolate the effect of φ(f T ) on front propagation from the influence of front initiation and stress state we leave the loading zone on the left unmodified (it is the same in both insets) so that the restarted event begins like the original; we also let the modified φ(f T ) have the same mean value as the original φ(f T ) for all blocks. Thus, the stress state is the same and the only change is in the width of φ. (A) Uniform and bell-shaped initial distributions of various widths. The bell-shaped distributions are polynomials with roots at ±a and the functional form φ(ξ) = 5/(4a)(1 + 3|ξ/a|)(1 -|ξ/a|) 3 , ξ ∈ [-a, a]. (B) Friction to normal force ratio vs block displacement for a block having the spring force distributions in A (corresponding colours). Spring relaxation during slip is excluded from the calculation, a valid assumption when the passage of the rupture front is quick compared to the mean slipping time t R . Markers are located at maxima, which define τ max /p. (C) Effective static friction coefficient vs σ, the standard deviation of φ. Each marker takes its abscissa from the data in A and its ordinate from B. Drawn and dashed lines connect markers corresponding to uniform and bell-shaped φ, respectively. They are used in Fig. 4C as reference lines for simulation data. For clarity, only a few of the φ that were used to determine the lines are shown on this figure. The form μeff s is analogous to τ0 introduced in Supporting Figure S3, and represents the ratio of (i) the dynamic stress drop from the effective static friction level to the stress obtained during sliding and (ii) the maximum value this stress drop can attain. All data in this figure comes from simulations with a model in which the slipping force f slip decreases linearly with time. This microscopic friction law modifies the slow slip mechanism with respect to the reference microscopic model. Until depinning, the springs have the same behaviour in both models. Upon entering the slipping state, the force from each spring on the block immediately drops to the level f 0 slip = µ d f N,ij , as in the reference model. Then the slipping force decays linearly in time with a decay rate depending on the slipping time and the force level at repinning f new ; that is, f slip = f 0 slip -(f 0 slip -f new )t s /t R , with t s measured from the time the spring entered the slipping state. The slipping time distribution is the same as in the reference model. The decay in the friction coefficient enables a slow slip motion after the initial rapid slip, which allows this model to sustain slow fronts for a relatively wide range of initial spring stretching configurations and stress states, as seen in B.

The parameters used with this model are the same as in the reference model (Supporting Table S1), except for k ij = 54.1 GN/m 2 f N,ij , µ d = 0.28 and the new parameter µ new = 0.7µ d ≈ 0.2 used to determine

f new = µ new f N,ij .

Supporting Equations

Estimate of v slow slip from model parameters

Here we provide the arguments behind the slow slip speed estimate used in the inset of Fig. 3C.

The mechanism for slow slip in our model is the force drop when slipping junctions relax and repin at zero force. To determine the slip speed (block motion per unit time) associated with the relaxation of junctions, we identify the net slip caused by this change and the time over which the change happens. To do this, we assume that each time a junction relaxes, the block moves forward just enough to return to the force it was bearing just before relaxation. Thus, the force reduction -f slip in the junction is matched by an equivalent net force increase in the junctions that are already pinned and in the coupling to the neighbouring blocks. The effective stiffness of these interface and bulk terms depend on the fraction of junctions that are pinned and on the motion of the neighbouring blocks, respectively.

A careful look at the motion of blocks during their slow slipping regime proved that neighbouring blocks move essentially at the same slow slip speed. This means that the force changes on a block due to relative motion with respect to its neighbours remain negligible. We can therefore assume in this calculation that the only contribution to force restoring after junction relaxations is due to the pinned interfacial junctions.

We now assume that all junctions start in the slipping state. After the first junction relaxes and repins, the effective stiffness of the interface is just the stiffness of this single junction, k ij . The force drop -f slip must be compensated by stretching this (now pinned) junction by moving the block a distance δ 1 = f slip /k ij . For the second junction the force drop is again -f slip , but the stiffness of the interface has increased to 2k ij and the required slip is δ 2 = f slip /(2k ij ). If no junction would reach its breaking threshold f thres during the whole relaxation process, then we would find x slow slip = Ns j=1 δ j = Ns j=1 f slip /(jk ij ). Since the k ij are independent of j this is just δ = f slip /k ij Ns n=1 1/n with n a dummy index. For N s = 100 used in the model, the sum evaluates to 5.2. However, for f thres /f slip = 40/17 used in the model, some junctions do break again before all slipping junctions relax. The force drop associated with the transition from pinned to slipping state is f slip -f thres , and acts in the same way as the force drop when junctions leave the slipping state. Taking this into account and evaluating the return to the pinned state more carefully, we find the net block slip to be 8.5f slip /k ij .

With the slipping time standard deviation δt R = 0.3 tR used in the model, the time for all the junctions to return to the pinned state is found to be close to 2 t R .

Combining these slip distance and slip time values, and defining τ slip = N s f slip and k i = N s k ij , we use

v slow slip, estimate = 4.2 τ slip k i t R . (1)
The inset of Fig. 3C plots the slow slip speed measured in the simulations against this estimate. The markers indicate when we have varied τ slip ( ), k i ( ) and t R ( ). • uses our reference parameters. The blue data is based on restarting event II. Because τ slip and k i enter in the elastic state of the slider, only t R could be varied for these simulations. The green data is for prepared homogeneous interfaces, where parameters can be varied freely. 1When slow slip speed is measured in the simulations based on event II (blue data in inset of Fig. 3C), some deviation from v slow slip, estimate are observed. There are several possible reasons for such a deviation. First, the assumption of co-moving neighbours is only approximately correct (near the front tip the neighbours to the right are stuck until the rupture front passes them). The actual During this time the front tip has moved the distance dx from one block to the next, and the front speed is

v c,slow = dx dt = v slow slip dx ∆u(x 0 ) (8) = v slow slip k i τ thres -τ 0 dx f [x 0 ] -f [x 0 + dx] . (9)
Here we recognize an approximation to the spatial derivative of the unknown function f (), evaluated at x = x 0 . We will use the chain rule to separate the non-dimensional and dimensional parts of this derivative, and therefore we define X(x) = (x -x 0 )/l 0 so that f

We arrive at

v c,slow = v slow slip k i τ thres -τ 0 l 0 -df dX x=x 0 . ( 11
)
This argument provides a rationale for the linear relationship observed in Fig. 3D, but with the function f () unknown we are not able to predict the value of the coefficient of proportionality. From the shear force profiles we estimate the decay length l 0 = 7 mm, a value shared between simulations because we keep k/k i the same, and rescale v slow slip with k i l 0 /(τ thres -τ 0 ). Note that in the model, τ thres = N s f thres = µ s p, with µ s the threshold force coefficient and p the normal force on the block, which means that the normal force enters in the scaling.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Sketch and behavior of the multiscale model. (A) Slider and external loading conditions. (B) Spring-block network modeling elastodynamics. (C) Surface springs modeling friction on a block. (D) Macroscopic loading curve, the ratio F T /F N of driving shear force to total normal force. (E) Mesoscopic loading curve, the ratio τ/p of shear to normal stress on a block. (F) Microscopic friction model for the spring loading curve, the ratio f T /f N of friction to normal force for one spring (f N = p/N s ). (G) Steady-state friction force on a block vs. sliding velocity.

FigFig. 3 .

 3 Fig.1D). ly event tion for e springs sured as B) using

Fig. 4 .

 4 Fig. 4. Role of microscopic force distributions in rupture dynamics. (A) Event II shown as in Fig. 2A. (B) Event II restarted at 1.054 s with wider ϕ(f T ) in the region x > 0.25L (arrow) results in a fast-only event. (C)μ eff s = ððτ max =pÞðf slip =f N ÞÞ= ððf thres =f N Þðf slip =f N ÞÞ vs. σ = σ=ðf thresf slip Þfor all blocks in events I (green) and II (red), and for prepared events (black) with homogeneous τ 0 = ððτ 0 =pÞðf slip =f N ÞÞ=ððf thres =f N Þðf slip =f N ÞÞ. •, slow front; ×, fast front; ○, loading zone (x < 2 cm). Solid (dashed) line is for uniform (bell-shaped) ϕ(f T ) (Fig.S4). (D) Observed front type for prepared interfaces vs. both τ 0 and σ. Dots represents data points.

A B 5 .

 5 Dynamic vs. quasistatic slow fronts. (Upper and Lower Left) Macroscopic loading curves. (Upper and Right) Fraction of pinned shown as in Fig. 2A. (A) Event I restarted. (B) Data from a simulation where driving is applied at h = 0. Gray scale represents the original data with no change to V. Color represents the data in which the driving speed is set to V = 0 during front propagation. After setting V = 0, we only show color when at least 5% of a block's springs are slipping.

B

  Supporting FigureS2: A fast-slow event arising spontaneously in a simulation using the simplified friction law with a single interface spring per block. The parameter v slow slip spring = 1.5 mm/s (Supporting Methods). (A) Spatio-temporal plot of the instantaneous friction to normal force ratio. (B) Rupture front speed v c . Block rupture is defined to occur when the interface spring depins. Front speed is measured as in Fig.2.Position along interface, x (m)Modified φSupporting FigureS4: Microscopic force distributions significantly affect rupture fronts. The data in Fig.4A,B repeated with spring force distributions shown in detail. (A) Event II shown as in Fig.2. (B) Simulation behind A restarted at 1.054 s with a wider distribution φ(f T ) of shear forces results in a fast-only event. Insets: For each block along the interface, a color coded histogram of φ(f T ) at 1.054 s. The vertical axis shows the force level in individual springs, which extends up to f thres . The level f thres is shown by the upper white line; it is different for each block because it varies with normal force. The lower white line corresponds to f T = 0 N. The color denotes the fraction of springs found at each value of f T using an arbitrary bin width. This means that (apart from normalization) each vertical slice in the insets shows the same type of data as Supporting FigureS5A. Offset data: Fraction of slipping springs at 1.054 s.

C

  Supporting Figure S5: Dependence of effective static friction coefficient µ eff s = τ max /p on the distribution φ of forces in the interface springs.

  Results from an alternative model with a different individual behaviour of the interface springs, showing the robustness of our results against changes in the interfacial law. (A) A spontaneously arising fast-slow-fast event, analogous to that in Fig. 2A. (B) Observed front type for prepared interfaces, analogous to Fig. 4D.

  [x] = f (X(x)). With this notation,

Table S1 :

 S1 Model parameters. Parameters above the horizontal line are used in the same way as in Ref.[START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF].

	Supporting Name	Symbol	Value
	Slider length (x)	L	140 mm
	Slider height (z)	H	75 mm
	Slider width (y)	B	6 mm
	Number of blocks	N x	57
		N z	31
	Slider mass	M	75.6 g
	Block mass	m	M/(N x N z )
	Young's modulus	E	3 GPa
	Bulk spring modulus	k	3BE/4
	Bulk spring length Damping coefficient	l η	L/(N x -1) = H/(N z -1) √ 0.1km
	Normal load	F N	1920 N
	Elastic foundation modulus	k f	k/2
	Driving spring modulus	K	4 MN/m
	Driving height	h	5 mm
	Driving speed	V	0.4 mm/s
	Threshold force coefficient	µ s = f thres /f N 0.4
	Slipping force coefficient	µ d = f slip /f N	0.17
	Number of interface springs per block N s	100
	Interface spring stiffness	k ij	39.2 GN/m 2 f N,ij
	Slipping time mean	tR	2 ms
	Slipping time standard deviation	δt R	0.6 ms
	Triggering region width	x trigger	22.5 mm
	Triggering region prestress	τtrigger	0.3
	Time step duration	∆t	2 • 10 -7 s
	Extra damping coefficient	α	η/40

When f slip is varied, f thres is varied proportionally in order to keep the fast-slow-fast nature of the rupture front. For the same reason, when ki is varied the prestress in the triggering region is also varied slightly.

Supporting information for

Slow slip and the transition from fast to slow fronts in the rupture of frictional interfaces

Supporting Figure S3: Spatial distributions of prestress τ 0 /p for prepared states of the interface. In the front triggering region on the left, the prestress is τ0 = τ0 p -

3 for all prepared states used in Fig. 3D, 4C and 4D. In the front propagation region the prestress is homogeneous along the interface, at a value varied between prepared states, here τ0 = 0.2 (black) and τ0 = 0.3 (blue). After initialization, all springs in the triggering region are depinned simultaneously. Initiating the events in this way, rather than by driving the system until rupture is triggered, ensures that the force drop / energy release in the triggering region remains the same between simulations.

The non-dimensional form τ0 of the prestress represents the ratio between (i) the stress in excess of the stress obtained during sliding and (ii) the maximum dynamic stress drop that results from rupture. It is analogous to the so-called S classically used in seismology and to the form defined in Ref. [START_REF] Trømborg | Transition from static to kinetic friction: Insights from a 2D model[END_REF]. motion of the neighbours also depends on the stress state and the triggering of the event. Second, the assumption of force re-balancing every time a junction changes state is probably too strong.

Scaling of v c,slow with v slow slip

The data collapse in Fig. 3D is obtained by plotting the slow rupture speed v c,slow against the quantity v slow slip k i l 0 /(τ thres -τ 0 ), where v slow slip is the slow slip speed, k i is the stiffness of the connection between a block and the interface (a single spring in the simplified model and a parallel connection in the reference microscopic model), l 0 is a characteristic length, τ thres = N s f thres is the maximum shear strength of a block and τ 0 is the shear force in the propagation region before the event is started. In this section we provide a crude argument for this scaling.

When a region of initially homogeneous prestress is being stressed further by block motion on the left, the decaying shear force profile can be written on the form

where A is an amplitude and f () is a function that has magnitude 1 at x = x 0 and decays over a characteristic length l 0 that depends on the bulk to interfacial stiffness ratio k/k i . The function f () is known in 1D (Ref. [START_REF] Amundsen | 1D model of precursors to frictional stick-slip motion allowing for robust comparison with experiments[END_REF], equation ( 46)); in 2D it can be measured in an elastostatic model, but its exact form is not required for the present argument.

In a static situation τ (x) is balanced by the friction forces in the interfacial springs. Ignoring the width of the spring force distribution, the block at x 0 is at its static friction threshold when the force on it from its neighboring blocks is τ (x 0 ) = τ thres , which gives A = (τ thres -τ 0 ). The next block to the right, at position x = x 0 + dx, then has

Here we have used the short-hand notation f [x] = f (x -x 0 )/l 0 .

As the front tip moves from the block at x 0 to the block at x 0 + dx, the force on this block from its neighbours increases to τ (x 0 + dx) = τ thres . It will be useful to rewrite this as τ (x 0 + dx) = τ thres = τ (x 0 ) = Af [x 0 ] + τ 0 . The change in force on the block at x 0 + dx is

Assuming a corresponding change in the friction force allows us to relate the force change to a displacement of the block, namely

In the next step we will need the displacement of the block at x 0 during the same time interval. As the blocks are at closely spaced points in a deforming elastic medium we will assume ∆u(x 0 ) = ∆u(x 0 + dx)(1 + O(dx)) ≈ ∆u(x 0 + dx). Now we make the approximation that after breaking, the blocks move at a constant speed v slow slip . It follows that the time it takes from when the block at x 0 breaks and until when the block at x 0 + dx breaks is dt = ∆u(x 0 ) v slow slip .


On Modeling Ego-Motion Uncertainty for Moving Object

Detection from a Mobile Platform Dingfu Zhou 1,2 , Vincent Frémont 1,2 , Benjamin Quost 1,2 and Bihao Wang 1,2 Abstract-In this paper, we propose an effective approach for moving object detection based on modeling the ego-motion uncertainty and using a graph-cut based motion segmentation. First, the relative camera pose is estimated by minimizing the sum of reprojection errors and its covariance matrix is calculated using a first-order errors propagation method. Next, a motion likelihood for each pixel is obtained by propagating the uncertainty of the ego-motion to the Residual Image Motion Flow (RIMF). Finally, the motion likelihood and the depth gradient are used in a graph-cut based approach as region and boundary terms respectively, in order to obtain the moving objects segmentation. Experimental results on real-world data show that our approach can detect dynamic objects which move on the epipolar plane or that are partially occluded in complex urban traffic scenes.

I. INTRODUCTION AND RELATED WORK

Vision-based driver assistance system (DAS) is a complex and challenging task in urban traffic scenarios. In particular, moving object detection from dynamic scene analysis is essential for obstacle avoidance and path planning, and has numerous applications in autonomous and semi-autonomous driving. Indeed, being able to detect dynamic obstacles (vehicles, cyclists, or pedestrians) and to estimate their positions and motion tendency can increase the safety in both autonomous and semi-autonomous driving. Moving objects information also helps to improve the precision of Visionbased Simultaneous Localization and Mapping (VSLAM) and Structure-from-Motion (SfM) approaches which mainly rely on static environment assumptions [START_REF] Bailey | Simultaneous localization and mapping (SLAM): part II[END_REF].

Several vision-based motion detection approaches have been proposed over the last decade. Using one camera, approaches like background subtraction [START_REF] Sheikh | Background subtraction for freely moving cameras[END_REF], adaptive background model [START_REF] Faro | Adaptive background modeling integrated with luminosity sensors and occlusion processing for reliable vehicle detection[END_REF] or optical flow measurement [START_REF] Talukder | Real-time detection of moving objects from moving vehicles using dense stereo and optical flow[END_REF] can be used when the camera is static. More details about motion detection from a static camera can be found in [START_REF] Nascimento | Performance evaluation of object detection algorithms for video surveillance[END_REF]. The problem becomes much more complex when the camera and the surrounding objects move simultaneously. Indeed, the camera motion induces location changes of all the image pixels. Therefore, geometrical constraints are essential to distinguish between static and moving parts of the image. Two-view geometrical constraints (known as epipolar constraints) can be used to detect moving pixels. However it cannot detect objects moving on the epipolar plane (degenerate case). Other constraints, such as flow vector bound constraints [START_REF] Abhijit Kundu | Moving object detection by multi-view geometric techniques from a single camera mounted robot[END_REF] and multi-frame epipolar constraints [START_REF] Dey | Detection of independently moving objects in non-planar scenes via multi-frame monocular epipolar constraint[END_REF], have been used to detect the objects moving on the epipolar plane. Using two cameras, a dense or sparse disparity map can be calculated to reconstruct 3D information of the The authors are with 1 Université de Technologie de Compiègne (UTC), 2 CNRS Heudiasyc UMR 7253, France. E-mail: {dingfu.zhou, vincent.fremont, benjamin.quost, bihao.wang}@hds.utc.fr environment [START_REF] Hirschmuller | Accurate and efficient stereo processing by semi-global matching and mutual information[END_REF]. By combining disparity information with feature tracking or optical flow computation, the 3D scene flow can be reconstructed and used to detect moving objects [START_REF] Lenz | Sparse scene flow segmentation for moving object detection in urban environments[END_REF], [START_REF] Wedel | Stereo scene flow for 3D motion analysis[END_REF].

Following [START_REF] Talukder | Real-time detection of moving objects from moving vehicles using dense stereo and optical flow[END_REF], we propose to detect moving using the Residual Image Motion Flow (RIMF), which quantifies the difference between the measured optical flow and global image motion flow (pixel changes caused by camera motion). The global image motion flow between two consecutive frames can be determined by a function of the current scene depth and the relative camera pose. Regions with significant RIMF are detected as potential moving objects. In order to avoid a large number of false positives or misdetections, the noise in the RIMF estimation process should be considered. In [START_REF] Wedel | Stereo scene flow for 3D motion analysis[END_REF] and [START_REF] Romero | Stereo-based motion detection and tracking from a moving platform[END_REF], the uncertainties of the real optical flow and 3D scene flow have been modeled respectively to detect the moving objects. However, they just roughly modeled the uncertainty of the ego-motion information obtained from other sensors (GPS/IMU).

Unlike these methods, our approach is only based on two consecutive stereo images: no other sensor information is required. Furthermore we detail how the ego-motion uncertainty may be taken into account so as to improve the RIMF computation. Therefore, in this paper we propose a moving object detection approach based on [START_REF] Zhou | Moving objects detection and credal boosting based recognition in urban environments[END_REF] with two main contributions. First, a first-order error propagation framework is used to take into account the ego-motion imprecision which result from the uncertainty in feature extraction and matching. The uncertainty propagation strategy is also applied in the RIMF computation to build a motion likelihood for each possible image pixel (when its disparity value is available). This information quantifies the likelihood for a pixel to be moving or not. Then, a segmentation of the moving objects is performed using graph-cuts based motion segmentation on the motion likelihood and depth information, which helps to reduce the noise in the local optical flow estimation process.

Our paper is organized as follows: first, we present the overview of our approach in Section II. Section III introduces the key steps of our moving object detection algorithm, including ego-motion estimation and uncertainty computation, motion likelihood calculation and graph-cuts based motion segmentation. Simulation experiments to test the uncertainty estimation and real experiments results for moving object detection are presented in Section IV. Finally, the paper ends with conclusions and future work. II. SYSTEM OVERVIEW Fig. 1 outlines the main steps of our moving object detection system based on two consecutive stereo image pairs. The dense disparity map at frame t -1 and the optical flow (dense or sparse) between frame t -1 and t are estimated respectively. At the same time, the feature points tracked and matched between two consecutive stereo image pairs are used to estimate the camera relative pose and its uncertainty.

Then, the global motion flow (image changes caused by the camera ego-motion) is calculated for each image pixel using the camera motion and the depth information obtained from the disparity map. The RIMF, which is used to measure the difference between the measured optical flow and the global image motion flow, is used to distinguish between moving and non-moving pixels. In order to handle the noise involved in the RIMF calculation, the covariance matrix of RIMF is used to estimate the motion likelihood for each image pixel. It is computed using a first-order strategy to propagate the uncertainty from both disparity and egomotion estimation procedures to the RIMF estimation. Using the motion likelihood and disparity values, a graph-cut based motion segmentation approach is then applied to segment the image into regions of moving and non-moving pixels. Finally, the bounding boxes of the moving objects are obtained by combining region-growing and U-disparity map information.

III. MOVING OBJECT DETECTION FRAMEWORK

Two successive stereo frames are considered in the motion detection procedure. We assume that the stereo rig undergone an unconstrained motion (R, T) between the two successive frames, where R and T are respectively the rotational and translational components of the motion. The left image at time t -1 is considered as the reference image. A 3D point P at frames t -1 and t is noted as P

t-1 (X t-1 , Y t-1 , Z t-1 ) and P (X t , Y t , Z t ) respectively. The selected 2D points (x L t-1 , y L t-1 ), (x R t-1 , y R t-1 ), (x L t , y L t ), (x R t , y R t )
are matched and tracked [START_REF] Geiger | Stereoscan: Dense 3d reconstruction in real-time[END_REF] ( in the left and right images at time t -1 and t (a bucketing technique is used to ensure that features points are well spread in the whole image regions). Assuming the origin of 3D coordinate system is coincident with the left camera center, the 3D world points can be obtained as follow:

(X t-1 , Y t-1 , Z t-1 ) T = b d (x L t-1 -u 0 , y L t-1 -v 0 , f ) T (1)
where d = x L t-1 -x R t-1 is the disparity value for the scanline y = y L t-1 . The variables f, b and (u 0 , v 0 ) are the camera intrinsic parameters known as the focal length, the baseline and the principal point coordinates.

A. Ego-motion Estimation and Uncertainty Computation 1) Ego-motion Estimation : Given the points matching in four images for two consecutive frames, the relative pose of the camera can be estimated by minimizing the sum of the reprojection errors using non-linear minimization approaches. First, the feature points from the previous frame are reconstructed in 3D via triangulation and using the camera intrinsic parameters. Then these 3D points are reprojected into current image frames using the camera motion as below:

xi t = f (Θ, x i t-1 ) = P r L (K [R|T] P i t-1 ) P r R (K [R|T] P i t-1 ) (2)
where

xt,i = (x L t,i , ŷL t,i , xR t,i , ŷR t,i ) T and x t-1,i = (x L t-1,i , y L t-1,i , x R t-1,i , y R t-1,i )
T are the predicted and measured image points at image t and t -1 respectively. The vector Θ = (r x , r y , r z , T z , T y , T z ) T represents the six degrees of freedom of the relative pose. Let P r L and P r R be the image projections of the 3D world points into the left and right image (non-homogeneous coordinates). Let P i t-1 be the 3D point in the previous frame, which is calculated using Eq. ( 1). The parameter vector Θ can be estimated using the minimization of the following cost function which is built using the geometric distance of the predicted and measured image points in time t as:

F (Θ, x) = x t -xt 2 Σx t = x t -f (Θ, x t-1 ) 2 Σx t (3)
where . 2 Σ denotes the squared Mahalanobis distance according to the covariance matrix Σ. A Gaussian-Newton iterative optimization method is used to solve the optimization problem presented in Eq. (3).

2) Error Propagation: For vision systems, robust motion estimation should not only provide an estimate of the camera motion, but also an estimate of the uncertainty associated with this solution. Let the i th loop matched features be

(x L t-1,i , y L t-1,i ), (x R t-1,i , y R t-1,i ), (x L t,i , y L t,i ), (x R t,i , y R t,i ), i = 1, • • • , N in stereo
image pairs in t -1 and t frames. By stacking these features in vectors, new vectors can be defined as: x ∈ R 8N represents all the features and x t ∈ R 4N , x t-1 ∈ R 4N represent the features at time t and t -1 respectively. To be robust against outliers (mismatched features or features on moving objects), a RANSAC strategy is applied to estimate the relative pose between the two successive frames. Assuming that all the inliers used for the final minimization of Eq. ( 3) are good matched features with only additive Gaussian noise, the associated probability density function has the following form:

x ∼ N µ xt-1 µ xt , Σ xt-1 0 0 Σ xt (4)
where, µ t and Σ t (respectively, µ t-1 and Σ t-1 ) are the mean and the covariance of the image features at time t (resp. t-1).

Then the parameters accuracy only depends on the precision of the detected feature locations in the image plane.

In [START_REF] Clarke | Modelling uncertainty: A primer[END_REF] and [START_REF] Robert | Propagating covariance in computer vision[END_REF], the covariance matrix of the estimated parameters which considers the uncertainty of x t and x t-1 respectively, can be obtained with the following model:

Σ Θ = ( ∂g ∂Θ ) -1 ( ∂g ∂x ) T Σ X ∂g ∂x ( ∂g ∂Θ -1 ) T (5)
where, g(x, Θ) = ∂F (X,Θ)

∂Θ

, is the partial derivative of F (Θ, x) w.r.t each component of Θ. The matrix Σ x , which has been defined in Eq. ( 4), is the covariance matrix of the measured features at time t -1 and t.

In order to get the partial derivatives, ∂g ∂Θ and ∂g ∂x of g(Θ, x) .w.r.t each components of Θ and x, and since ∂g ∂Θ and ∂g ∂X are parts of the Hessian of F (Θ, x), we have:

∂g ∂Θ = ∂ 2 F ∂Θ 2 = 2( ∂f ∂Θ ) T Σ -1 x t ∂f ∂Θ -2 ∂ 2 f ∂Θ 2 Σ -1 x t (xt -f (Θ, xt-1)) (6)
∂g ∂x = ∂ 2 F ∂Θ∂x = ∂ 2 F ∂Θ∂xc ∂ 2 F ∂Θ∂xt-1 = -2Σ -1 xt ∂f ∂Θ ∂ 2 F ∂Θ∂xt-1 (7) ∂ 2 F ∂Θ∂xt-1 =2( ∂f ∂xt-1 ) T Σ -1 x t ∂f ∂Θ -2 ∂ 2 f ∂Θ∂xt-1 Σ -1 x t (xt-f (Θ, xt-1)) (8)
Equation. ( 6) and ( 8) are both a sum of two terms, a first part with only first-order derivatives and a second part which is the product of second-order derivatives and the residual part x t -f (Θ, x t-1 ). If the final solution of Eq. ( 3) has a zero residual or very small residual, one can remove the secondorder derivatives parts from Eq. ( 6)-( 8) to get a first-order estimate as below:

Σ Θ = ( ∂f ∂Θ T Σ -1 Xc ∂f ∂Θ ) -1 + A T Σ Xp A (9)
where,

A = 1 2 ∂f ∂Θ T Σ -1 Xc ∂f ∂Θ -1 ( ∂f ∂Θ ) T Σ -1 Xc ∂f ∂X p
Note that using only the first right-hand term in Eq. ( 9) is called as partially method (only the current frame feature noise has been considered) which has been described in [START_REF] Zisserman | Multiple view geometry in computer vision[END_REF]. However, due to the noise of the measured features in the previous and the current frames, the residuals of the final solution of Eq. (3) may not be zero. So a second-order error propagation model which considers residuals will be more suitable. A second-order error propagation result can be obtained by substituting Eq. (6,8) into Eq. ( 5).

B. Residual Image Motion Flow

Given a pixel position x t-1 , it is possible to compute the image position of pixel xt at time t using the following equation [START_REF] Zisserman | Multiple view geometry in computer vision[END_REF]:

xt = KRK -1 x t-1 + Kt z (10)
where x t-1 and xt are the normalized homogeneous image pixels at time t -1 and t respectively. The matrix K encapsulates the camera intrinsic parameters, z is the depth of the image pixel x t-1 . Normalizing Eq. ( 10) and substituting z = bf d into Eq. ( 10), the following equation can be obtained:

x t = r00f b(xt-1-u0)+r01f b(yt-1-v0)+r02f 2 b+df Tx r20b(xt-1-u0)+r21b(yt-1-v0)+r22f b+dTz + u 0 y t = r10f b(xt-1-u0)+r11f b(yt-1-v0)+r12f 2 b+df tTy r20b(xt-1-u0)+r21b(yt-1-v0)+r22f b+dTz + v 0 (11
) where r ij , i, j = 1, 2, 3 are the rotation matrix coefficient in line-column format. For image pixels from 3D static points, the global image motion flow caused by the camera motion only can be expressed as:

u x v y = x t -x t-1 y t -y t-1 (12)
However, for moving objects, the image motion is caused by both their motion and the camera displacement. Assuming that the real optical flow estimated from the image between time t and t -1 is (u x , v y ), it is possible to define the RIMF as:

--→ p (m) = (u x -u x , v y -v y ) (m) (13)

C. Motion Likelihood Estimation

Once the RIMF has been computed using Eq. ( 13), it can be used to separate the image pixels into moving or non-moving parts. Comparing the absolute RIMF difference | ---→ p (m) | to a fixed threshold does not lead to a satisfying result to differentiate moving pixels from static ones. Points with different 3D world locations have different image motions. Also, the estimation uncertainty, e.g. camera motion and pixel depth, have different influences on the image points. Ignoring these uncertainties could lead to a large number of false positives. In our case, the uncertainty in the RIMF is propagated from the image pixels noise to the final estimation using a first order Gaussian approximation.

As in Eq. ( 13), the RIMF is a function of the camera motion Θ, the pixel location (x t-1 , y t-1 ) in the previous frame, the depth of its corresponding 3D point d and the measured optical flow (u x , v y ). The uncertainty of the measured optical flow will not be considered in this work because it only affects the detection result locally. However, a linear approximation of the RIMF covariance can be calculated as:

Σ RIM F = JCJ T
where J represents the Jacobian matrix with respect to the camera motion Θ, the pixel position (x t-1 , y t-1 ) in the previous frame and the disparity value d in previous frame, and C is the covariance matrix of all the input variables: with Σ Θ the covariance matrix associated to the camera motion estimated in Sec.III-A and Σ 1 = diag σ 2 x , σ 2 y , σ 2 d . The variances σ x and σ y describe the image noise coming from features detection. As in [START_REF] Wedel | Stereo scene flow for 3D motion analysis[END_REF], the disparity uncertainty can be considered as an approximate standard Gaussian Distribution and the variance can be approximated by a linear function,

σ d (x, y) = σ 0 + γU d (x, y) (14)
where σ 0 and the γ are two constants parameters and U d (x, y) is the uncertainty on the disparity value at position (x, y). Here, the matching cost is used as a confidence measure of the disparity value (further details can be found in [START_REF] Hu | A quantitative evaluation of confidence measures for stereo vision[END_REF]).

Assuming a stationary world and Gaussian error propagation, static pixels can be expected to follow a Gaussian distribution with zero mean and covariance matrix Σ RIM F . The Mahalanobis distance associated to the RIMF can be calculated as follows:

µ p (m) = p T m Σ -1 RIM F p m (15)
Sinceµ 2 p (m) isχ 2 distributed, the motion likelihoodξ motion (x) of the RIMF of each pixel can be computed according to its µ p (m) value.

D. Graph-based Segmentation

Graph-cut (GC) is an energy minimization framework widely used in image segmentation. Further details may be found in [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in ND images[END_REF] [START_REF] Boykov | Graph cuts and efficient ND image segmentation[END_REF]. Usually independent moving objects are layered in the depth direction in real traffic scenes and may thus easily be distinguished from the neighboring background. Here, GC is applied to depth images in order to refine the results obtained from likelihood motion estimation. In particular, we detail how to build our cost function for segmentation, which is composed of two terms related to local and contextual segmentation:

E(L) = E r (L) + λE b (L) whereL = {l 1 , l 2 ,
• • • , l p }is a binary vector,l i ∈ {0, 1}is the label of the pixel (1 if moving, and 0 if static). Here, E r and E b are called the regional and boundary terms, respectively. The former term E r (L), ensures that pixels with a high motion likelihood are to be detected as moving; the latter, that adjacent pixels with similar depths should share the same label. The parameter λ is used as a trade-off between both terms. 1) Region Term: The motion likelihood of each pixel can be used directly to build the region term E r of the energy function.

E r = - Ω {L(x)ξ motion (x) + (1 -L(x))ξ static (x)}
where Ω represents all the image field and ξ static is a fixed prior likelihood to describe a point to be static. We set ξ static = 0.5 in our experiments.

2) Boundary Term : Usually, the depth maps can be used for object segmentation [START_REF] Hernã¡ndez-Vela | Graph cuts optimization for multi-limb human segmentation in depth maps[END_REF] because the depth of objects is significantly different from the background. Therefore, we propose to use the depth de(x) = bf D(x) (D is the disparity value) to measure the similarity between two pixels x i and x j as follows:

B(x i , x j ) = exp(- √ 2(|de(x i ) -de(x j )|))
Then, the boundary term can be expressed as below:

E b = Ω x∈N4(x) B(x i , x i )|L(x) -L(x)|
where N 4 (x) is the 4-neighborhood of a pixel x. In order to improve the segmentation efficiency, down-sampling is used in the motion likelihood estimation and moving objects segmentation steps, retaining one pixel out of four in both dimensions of the image.

E. Bounding Box Generation and Verification

After the segmentation step detailed above, bounding boxes should be generated for every moving object before performing tracking or recognition. Note that errors may come from partial detection (e.g., legs or arms of pedestrians) or redundancies (such as shadows). Object verification and region growing may be used to remove redundancies and to integrate parts detection using the dense disparity map. The U-and V-disparity maps [START_REF] Hu | UV-disparity: an efficient algorithm for stereovision based scene analysis[END_REF] are two variants of the disparity map that are widely used for road and obstacle detection. In the U-disparity map, an upright object will form a horizontal line because of the same disparity value. This information may be used to obtain the width of the bounding box. Then, region growing can be applied to get the height of the bounding box from the disparity map. According to [START_REF] Hoiem | Putting objects in perspective[END_REF], the real world height of the objects could be estimated as below:

h i = h c + (y i -y 0 )z cos θ f (16)
Here, h i and h c are the height of object i and camera respectively in real world coordinate; θ is the camera tilt angle and f is the camera focal length; z is the depth of the object to the camera; y 0 and y i are the horizon position and top of the objects in image coordinate (the origin of coordinates is assumed at left bottom). Assuming that moving objects are not higher than 2.5 meters, some obvious false positives may be filtered. For this purpose, the horizontal position is first computed using the V-disparity map. Then, the actual height of the object h i is calculated using Eq. ( 16). Finally we retain only the objects which height is between 0.5m and 2.5m.

IV. EXPERIMENTS

A. Motion Uncertainty Estimation 1) Monte Carlo Experiments : In our simulation experiments, both intrinsic and extrinsic parameters of the stereo rig are assumed to be known. The relative pose of the stereo cameras between two successive frames is fixed before generating the image features. We generate 3D points from a uniform distribution. Then they are projected into the four images using the appropriate projection matrices. We use bucketing techniques in order to ensure that features points are well spread in all the image regions.

A Monte-Carlo-like experiment is used to obtain an estimate of the covariance matrix as ground truth. At each time, the measured features are generated using Eq. ( 4) are used as inputs in Eq. ( 3) to obtain the optimal parameters in Θ . Covariance matrices of the Monte-Carlo method can be calculated from N independently estimated Θ. We set N = 500 in our experiments.

Simulation experiments have been conducted to compare the performance of the different approaches. Fig. 2 clearly 5) respectively, and partially method (blue) using in [START_REF] Zisserman | Multiple view geometry in computer vision[END_REF].

shows that the first and second-order methods proposed in Sec.III-A perform better than the classical partially method technique, the latter being slightly superior to the former; both perform almost as well as the Monte-Carlo approximation. For the sake of computational efficiency, we used first-order method in all our experiments. In order to validate our moving object detection approach, we used the video sequences from the KITTI dataset [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF] with a resolution of 1240×370 pixels and 10 frames per second. The video sequences are acquired from a SVS installed on the roof of a vehicle. Details about the setup can be found in [START_REF] Geiger | Are we ready for autonomous driving? the KITTI vision benchmark suite[END_REF]. Five video sequences acquired in different road situations were used to test our moving object detection technique. In the inner city sequence, the host vehicle was driven at a low speed (15 km/h) because of the complex road conditions. In the suburban road, the speed went up to 50 km/h. First, the stereo disparity map [START_REF] Geiger | Efficient largescale stereo matching[END_REF] and optical flow (dense [START_REF] Liu | Beyond pixels: exploring new representations and applications for motion analysis[END_REF] or sparse [START_REF] Geiger | Stereoscan: Dense 3d reconstruction in real-time[END_REF] ) are computed before the moving object detection step. At the same time, the relative pose of the camera between two consecutive frames and its covariance are estimated as mentioned in III-A, and we assumed that the covariance of the features in Eq. ( 4) is 0.5 pixel. To compute the variance of disparity in Eq .( 14), we set σ 0 = 0.25 and γ = 0.075 empirically. Fig. 3 shows the key steps of our moving object detection approach. Fig. 3(a) is the motion likelihood image, the stationary and moving parts are respectively displayed in green and red. Fig. 3(b) and 3(c) show two detections results based on a fixed threshold and graph-based method. From the results, we can see that our graph-cut-based approach performs better to detect moving objects than when using a fixed threshold. Note that the verification in the bounding box generation step allowed to avoid a false positive detection. Fig. 4 shows the detection results in a campus sequence. During this sequence, the camera turned from left to right at a high speed. Our algorithm proved to be efficient in this context too, allowing to detect the cyclists behind the trees to the left. We also tested our algorithm on a suburban road; the detection results are displayed in Fig. 5. In this sequence, both the camera and the object vehicle move at urban speed (about 50km/h). In this sequence, we used sparse instead of dense optical flow because of the high changes in the images between two successive frames. The opposite driving vehicles were detected at a range of 40m, which remains sufficient for an appropriate reaction of the driver. The white car moving in front of the camera was also properly detected even if it moves in the same direction that Table I: Performance of our moving object detection and segmentation approach the ego vehicle. Fig. 6 displays the results obtained on three inner city sequences. In crowded streets, the host vehicle moved slowly, which makes detecting moving objects easier, even when they move on the epipolar plane. Note that the algorithm also detected partially occluded objects because of the use of a dense approach. Despite the effectiveness of our algorithm, some false negative and false positive detections happen in the real image sequences. In Fig. 4, a cyclist (red elliptical box) has not been detected because the related 3D information cannot be reconstructed in this frame (the cyclist can not be seen in right camera). False positives also appeared, such as in Fig. 6 (a) with red bounding box, due to reflections on windows in the scene. Table I describes our detection results in one inner city sequence with 153 frames (Fig. 6 (b)). The ground truth of the object's locations have been included in the dataset for this sequence. Here, only the moving objects whose distance is less than 30m are considered. The detection result is only based on two adjacent frames. True moving represents the number of detected moving object bounding boxes in the whole sequence. False moving means static objects that have been detected as moving and false static are the moving objects that are not detected. The true static objects are not taken into account because our algorithm focuses on moving objects only. In this sequence, our algorithm obtains a precision of 79.5% along with a recall of 86.7%.

B. Moving Object Detection

All the experiments have been realized on a standard laptop (Intel Core i7) with Matlab R2013a processing environment. When the dense optical flow is used, the total average computational time is about 30 seconds for each frame. The optical flow calculation step takes about 15 seconds. Around 4.5 seconds is spent on the motion likelihood computation, 5 seconds on the graph-cut based segmentation and 5 seconds on the bounding boxes generation. Computing ego-motion and estimating the uncertainty only takes about 0.2 seconds. Although our Matlab implementation is not real-time, it compares favorably with respect to [START_REF] Kumar Namdev | Motion segmentation of multiple objects from a freely moving monocular camera[END_REF] (7 minutes per frame) and further accelerations could be achieved by C/C++ implementation with parallel/GPU computing.

V. CONCLUSIONS AND FUTURE WORK

In this paper, a novel approach has been proposed to detect moving objects from two consecutive stereo frames by modeling the ego-motion uncertainty and using a graphcut based motion segmentation. The ego-motion uncertainty estimated through a first-order error propagation model is used to obtain the motion likelihood for each image pixel. Pixels with a high motion likelihood and a similar depth are detected as a moving object based on a graph-cut motion segmentation approach. Additionally, a fast recognition of moving objects becomes possible based on our segmentation results. Detection results in several different real video sequences show that our proposed algorithm is highly robust with respect to global (camera motion) and local (optical flow) noise. The ego-motion error has been considered using its covariance matrix and uncertainties in optical flow can be eliminated by graph-cut segmentation procedure. Furthermore, our approach works with all image pixels and arbitrarily moving objects (including partially occluded ) can be detected.

Future work will be firstly to consider a robust multiple objects tracking, using for example a PHD Filter, to obtain a stable detection results by reducing false positive detections. Further, fusing detection results coming from other sensors (lidar or radar) will also be tested to improved the detection results. Furthermore, categories information like pedestrian, car or others about the moving objects can be used on each bounding box to focus on moving pedestrians only.

Figure 1 :

 1 Figure 1: Moving objects segmentation framework overview.

Figure 2 :

 2 Figure 2: Motion Uncertainty estimates. Covariances between (a) ty and rx and (b) tx and ty, estimated using Monte-Carlo (green), first-and second-order techniques (black and red) using Eq. (5) respectively, and partially method (blue) using in[START_REF] Zisserman | Multiple view geometry in computer vision[END_REF].

  (a) Motion likelihood image (b) Segmentation result using a fixed threshold (c) (c) Graph-based segmentation result using depth (d) Bounding boxes generation and verification

Figure 3 :

 3 Figure 3: Moving Object Detection Flow . (a) motion likelihood (red: moving, green: static). (b) Segmentation using a fixed threshold ξ motion = 0.75. (c) Graph-based segmentation using depth information, We choose λ = 1 and ξ static = 0.5 for all our experiments. (d) Final result after bounding box generation and verification.

Figure 4 :

 4 Figure 4: Detection results on a campus sequence.

Figure 5 :

 5 Figure 5: Detection results on a suburban road.

Figure 6 :

 6 Figure 6: Detection results on three different inner city sequences. (a) crowded street (b) relatively simple environment (c) intersection