Guillaume Rongier

Pauline Collon

Philippe Renard

Stochastic simulation of channelized sedimentary bodies using a constrained L-system

Keywords:

Simulating realistic sedimentary bodies while conditioning all the available data is a major topic of research. We present a new method to simulate the channel morphologies resulting from the deposition processes. It relies on a formal grammar system, the Lindenmayer system, or L-system. The L-system puts together channel segments based on user-defined rules and parameters. The succession of segments is then interpreted to generate non-rational uniform B-splines representing straight to meandering channels. Constraints attract or repulse the channel from the data during the channel development. They enable to condition various data types, from well data to probability cubes or a confinement. The application to a synthetic case highlights the method's ability to manage various data while preserving at best the channel morphology.

I

The presence of channelized sedimentary bodies constitutes a key structuring element for the connectivity of a reservoir. These bodies are composed of heterogeneous deposits that can contain fluids or act as flow barriers. Such barriers may compartmentalize the reservoir and make its exploitation harder, especially when combined with sealing faults [e.g., [START_REF] Gainski | Turbidite reservoir compartmentalization and well targeting with 4D seismic and production data: Schiehallion Field, UK[END_REF]. Channel modeling can help to study and anticipate more precisely the flow behavior. It requires a balance between geological concepts and data conditioning, which remains an issue in stochastic simulation.

Stochastic simulation methods are usually split in two categories: cell-based and object-based methods. Cell-based methods [e.g., [START_REF] Deutsch | GSLIB: Geostatistical Software Library and User's Guide[END_REF][START_REF] Galli | The Pros and Cons of the Truncated Gaussian Method[END_REF][START_REF] Strebelle | Conditional simulation of complex geological structures using multiple-point statistics[END_REF] attribute a sedimentary body type -e.g., channel -to each cell within a grid, based on a prior model. It makes well data conditioning easy, but the channel continuity is rarely preserved. The resulting connectivity differs from the prior model or from object-based methods [START_REF] Rongier | Comparing connected structures in ensemble of random fields[END_REF]. Object-based methods [e.g., Viseur, 2001[START_REF] Deutsch | FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems[END_REF][START_REF] Pyrcz | ALLUVSIM: A program for eventbased stochastic modeling of fluvial depositional systems[END_REF] rely on a geometrical representation of the channels, with parameters such as their width or their wavelength. They preserve the channel continuity. However, data conditioning is difficult, because of the elongated shape of channels and the poor flexibility of their representations.

In numerical biology, tree and root simulation also relies on object-based approaches [e.g., [START_REF] Prusinkiewicz | The Algorithmic Beauty of Plants[END_REF][START_REF] Leitner | A dynamic root system growth model based on L-Systems[END_REF][START_REF] Longay | Treesketch: interactive procedural modeling of trees on a tablet[END_REF]. A formal grammar -the Lindenmayer system or L-system [START_REF] Lindenmayer | Mathematical models for cellular interactions in development I. Filaments with one-sided inputs[END_REF]] -mimics tree and root growth to simulate the related objects. An interesting development of such methods introduces the environment influence to improve the simulation realism: influence of gravity, influence of the sun light distribution, influence of other trees, etc. [e.g., [START_REF] Mvech | Visual models of plants interacting with their environment[END_REF][START_REF] Streit | Modelling Plant Variation Through Growth[END_REF][START_REF] Taylor-Hell | Incorporating Biomechanics into Architectural Tree Models[END_REF]. Despite this ability to integrate various data, very few works use formal grammars to simulate channels. [START_REF] Hill | Describing and generating facies models for reservoir characterisation: 2D map view[END_REF] rely on an analog model to define some rules for channel stochastic simulation, with a conditioning limited to well data or channel segments interpreted on seismic data.

In this paper, we present a new method to simulate channels based on L-systems. L-system rules control the channel morphology to simulate straight to sinuous channels (section 1). During the simulation, constraints influence the channel growth to condition various data (section 2). The purpose of this research is to facilitate data conditioning as compared to other object-based approaches while keeping their main addedvalue: the preservation of the channel shape. A synthetic case study applies the method to the simulation of submarine channels within an incised valley (section 3). It leads to a discussion on the channel simulation with L-system and constraints (section 4).

L--

This object-based method relies on a L-system to give form to a parameterized channel geometry.

Channel object definition and simulation principle

A channel object is built upon Non-Uniform Rational B-Splines (NURBS) -a generalization of Bézier curves [START_REF] Piegl | The NURBS Book[END_REF] -as defined by [START_REF] Ruiu | Modeling channel forms and related sedimentary objects using a boundary representation based on non-uniform rational B-splines[END_REF]. In this model, a channel appears as a discretized object composed of channel segments separated by transversal channel sections (figure 1).

A section is at the end of its respective segment. An orientation and a distance from the previous section characterize a given section. A L-system builds the succession of sections and determines their location. A sequential Gaussian simulation [Deutsch and Channel discretization for L-system simulation. The red dots locate the channel sections obtained from the interpretation. l s is the distance between two successive sections, i.e., the channel segment length. α is the angle between two successive sections, i.e., the change of orientation between two successive channel segments. The other parameters characterize a channel bend, with l B the bend length, λ the bend half-wavelength, ∆ the bend amplitude, and r = 1/c the radius of curvature -radius of a circle fitted at the bend apex -with c the curvature. [START_REF] Deutsch | GSLIB: Geostatistical Software Library and User's Guide[END_REF] simulates a width and thickness for each channel section, with a possible influence of the channel curvature. The section locations, widths and thicknesses are the foundation for the NURBS representation. This method is designed for a classical modeling workflow, which consists in building a geological grid from a structural model, and simulating the channels inside the parametric space of this grid [START_REF] Dubrule | Reservoir Geology Using 3-D Modelling Tools[END_REF].

Brief introduction to L-systems

The Lindenmayer system is a formal grammar designed by [START_REF] Lindenmayer | Mathematical models for cellular interactions in development I. Filaments with one-sided inputs[END_REF] to simulate the development of filamentous organisms. Since then, it has been expanded to simulate the development of plants [e.g., [START_REF] Mvech | Visual models of plants interacting with their environment[END_REF][START_REF] Prusinkiewicz | The use of positional information in the modeling of plants[END_REF][START_REF] Leitner | A dynamic root system growth model based on L-Systems[END_REF]. A L-system aims at rewriting an initial string, the axiom, with production rules, all composed of letters from a predefined alphabet. The rules include a set of letters to replace, the predecessor (pr ed), and a set of replacement letters, the successor (succ):

l c < pr ed(par am) > r c : cond p -→ succ(par am )
The application of a rule may depend on the letters before and after the predecessor, called respectively the left and right context (lc and r c), on a probability (p), and on parameters (par am and par am ) and conditions (cond). The following example contains an axiom ω and three production rules p 1 ,

Letters(parameters) Interpretation

I

First letter of a L-system string that contains the initial position

ω : b(0) p 1 : a(h) > b(h) 0.75 -→ b(h + 1)
p 2 : a(h) > b(h) 0.25 -→ a(h)a(h) p 3 : b(h) : h < 1 -→ a(h)b(h)
The rules may rewrite the axiom as follows:

µ 0 : b(0) (ω) µ 1 : a(0)b(0) (p 3 ) µ 2 : b(1)a(0)b(0) (p 1 & p 3 ) µ 3 : b(1)a(0)a(0)a(0)b(0) (p 2 & p 3 )
The resulting string constitutes a sequence of commands, which control an interpreter called a turtle [START_REF] Prusinkiewicz | Graphical applications of L-systems[END_REF]. The turtle's state is represented by a position vector P and an orthonormal coordinate system centered on this position, with:

• Ĥ the forward direction.

• L the left direction.

• Û the up direction. This state is built from an initial position and orientation. Then, the letters update it and draw object elements, such as segments, along the way to progressively build the object. More details about L-systems can be found in [START_REF] Prusinkiewicz | The Algorithmic Beauty of Plants[END_REF].

Alphabet for channel simulation

A channel is decomposed into bends separated by inflection points (figure 1). A bend is a succession of channel segments whose orientations change along the same direction. Seven letters constitute the alphabet to model this succession of segments (table 1).

A channel divides into two branches that grow in opposite directions. Brackets [ and ] symbolize this branching structure, with the letters of the first branch in between. The letter C builds the channel object, while the letters + andcontrol the channel sinuosity. A channel segment is symbolized by a ±C pair. Thus, a channel is a succession of ±C letters, with ± being + orbut remaining the same along the bend. The letter I has no influence on the string rewriting. It only states the initial position during the geometrical interpretation. The letter T has no influence on the geometrical interpretation. It only ensures the channel growth during the rewriting.

Production rule definition

The rules aim at building bends from channel segments. This includes converting the input parameters that describe the channel morphology into L-system parameters.

From input to L-system parameters

Simulating the channel morphology calls for:

• The global direction D along which the channel must be oriented.

• The default length l d for the channel segments to define the discretization resolution.

• The bend half-wavelength λ and the bend amplitude ∆ (figure 1) to define the bend morphology. They are drawn from statistical distributions for each new bend.

• The deviation angle δ to perturb the bend morphology, because a half-wavelength and an amplitude are not enough to characterize a bend [START_REF] Howard | Multivariate characterization of meandering[END_REF]. It is drawn from a statistical distribution for each new channel segment.

The bend length and the curvature can replace the bend half-wavelength, the bend amplitude and the deviation angle [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF] 1 . The purpose is then to retrieve the two Lsystem parameters l s and α.

It starts by determining the curvature c with the intersecting chord theorem and the law of cosines:

c = 8∆ 4∆ 2 + λ 2 (1)
c provides the bend length l B , i.e., the curvilinear distance between two inflection points (figure 1):

l B = 2 arccos(1 -∆c) c (2)
l B and l d determine the number n s of bend segments:

n s = l B /l d (3)
where x rounds x upward. The segment length l s becomes:

l s = l B /n s (4)
Finally, the angle α between two successive segments i -1 and i comes from the curvature and the lengths of both segments, perturbed by the deviation angle:

α = c(l s,i-1 + l s,i ) 2 + δ (5)

Rules for channel simulation

The channel initial position is randomly drawn inside the geological grid, possibly influenced by a grid property. The initiation rules simulate a single bend separated into the two branches (figure 2). The numbers of segments for each branch, n s 1 and n s 2 , are computed by:

n s 1 ∼ (1, n s -1) n s 2 = n s -n s 1 (6)
with n s the number of bend segments and a uniform discrete distribution define by its minimal and maximal values. n s 1 is randomly drawn from this distribution. Instead of using equation 5, the angles α 1 and α 2 of the first segment of each branch are modified into to initially orient the channel along the global direction:

α 1 = α 0 × ( n s 2 -n s 1 + 1) α 2 = α 0 × ( n s 2 -n s 2 ) α 0 = cl s + δ (7)
where δ is the deviation angle, c is the curvature (eq. 1), and l s is the segment length (eq. 4). The initiation rules are stochastic, with two rules for the two possible bend orientations. A letter T at each branch extremity enables further growth.

The channel growth is ensured by the development rules, which tie bends one after the other. These rules replace a T by a new bend depending on the left context, and add a T at the end to pursue the channel development (figure 2). This enables to change the orientation from one bend to the other. The process stops as soon as the two branches go outside the geological grid. The simulation can also stop at a given channel length. The complete rules for channel simulation are available in [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF]. These rules are predefined, but nothing prevents from modifying them.

Extending the influence of the global direction

Only the first segment of each branch aligns the channel along the global direction (eq. 7). When the L-system parameters do not vary much between bends, the channel keeps following the channel direction. When the L-system alternates straight and sinuous bends, the channel diverts from the global direction and often self-intersects. Those potential inconsistencies result from the L-system interpretation principle: a new segment only knows its immediate previous neighbor.

Extending the influence of the global direction to the whole system better controls the channel orientation and limits selfintersections. It occurs after the orientation change of the letters + and -. The global direction D reorients the turtle's head Ĥ:

Ĥ = K K K = ε h Ĥ + ε d o D (8)
with Ĥ the new turtle's head, ε h a weight on the L-system, ε d a weight on the global direction and o the branch orientation, equal to 1 for one branch and -1 for the opposite branch. ε h and ε d are user-defined, and control the influence of the Lsystem and the global direction.

C

In biology, constraints make trees interact with their environment [e.g., [START_REF] Taylor-Hell | Incorporating Biomechanics into Architectural Tree Models[END_REF][START_REF] Longay | Treesketch: interactive procedural modeling of trees on a tablet[END_REF]. Here, constraints condition the data.

Axiom IT Initiation Rules T 0.5 -→ +(α 1 )C{+C} ns 1 -1 T + (180) -(α 2 )C{-C} ns 2 T T 0.5 -→ -(α 1 )C{-C} ns 1 -1 T + (180) + (α 2 )C{+C} ns 2 T
Final string Iteration 1

I -(α 1 )C-C-C T + (180)+(α 2 )C+C T Interpretation Development Rules -C < T -→ {+C} ns T +C < T -→ {-C} ns T
Final string Iteration 2

I -(α 1 )C-C-C+C+C+C+C T + (180)+(α 2 )C+C-C-C-C-C T Interpretation Final string Iteration 3 I -(α 1 )C-C-C+C+C+C+C-C-C-C-C T + (180)+(α 2 )C+C-C-C-C-C+C+C+C+C T Interpretation Figure 2
Simplified illustration of channel simulation with a L-system. The red point is the initiation point. The first iteration stochastically initiates the two channel branches. The following iterations grow the channel. Unless otherwise specified, + andstand for +(α) and -(α) (see eq. 5), and C for C(l s ) (see eq. 4). α 1 and α 2 come from equation 7.

Formalization

An environmentally-sensitive process adds the constraints to the L-system. After each rewriting step, an interpretation defines the spatial relationship between the channel and its potential constraining elements. The constraints are separated into:

• Relative constraints, whose magnitude depends on the distance to their constraining element.

• Global constraints, whose magnitude is independent from a constraining element.

A constraint is also attractive or repulsive (figure 3).

Building constraints means choosing one or several constraining elements within a perception area (figure 4) and determining the constraint vectors. These vectors indicate where the system should go under the constraint influence. They decompose into three elements:

• The constraint direction Λ, defined from the vector V between the turtle's head and a point on the constraining element, either its center or the closest location to the turtle's head. (figure 3):

Λ =    V V if the constraint is attractive D-V D-V if the constraint is repulsive (9)
with D the global channel direction.

• The constraint magnitude ρ. It is equal to one for a global constraint, and depends on the distance to the constraining element for a relative constraint (figure 5).

• A user-defined weight ε to adjust the relative importance of each constraint. Constraint addition to a L-system, with Ĥ the current turtle's head, Ĥ the new turtle's head influenced by the constraint, D the channel global direction, Λatt the direction for an attractive constraint, deduced from V , Λrep the direction for a repulsive constraint, deduced from -V and D, and V the vector between the turtle's head and the constraining element.

Examples of constraining data

Wells provide local and often precise data about facies and sedimentary bodies. We consider two types of well data (table 2):

• Channel data comes from channelized bodies. A channel is attracted by one channel datum at a time until conditioning, i.e., the datum ends up in a channel.

• Inter-channel data correspond to other sedimentary bodies around channels, such as levees. A channel is repulsed by inter-channel data, i.e., these data must never end up in a channel.

Seismic data provide an overview of the underground in three-dimensions. We consider two constraints from seismic data (table 2):

• With enough resolution, seismic data can help interpret the channel confinement [e.g., [START_REF] Allen | Free Meandering Channels and Lateral Deposits[END_REF][START_REF] Janocko | The diversity of deepwater sinuous channel belts and slope valley-fill complexes[END_REF], a major source of uncertainty [START_REF] Larue | Why is reservoir architecture an insignificant uncertainty in many appraisal and development studies of clastic channelized reservoirs?[END_REF]. Here a geological grid materializes the confinement and constraints prevent the channels from going out (table 2).

• With a low resolution, a sand probability cube can be extracted [e.g., [START_REF] Strebelle | Modeling of a deepwater turbidite reservoir conditional to seismic data using principal component analysis and multiple-point geostatistics[END_REF]. As channels often have a sandy filling, the high sand probabilities attract the channels (table 2). The mean of many realizations (the E-map) should come close to the probability cube, without forgetting that simulating channels is not simulating sand distributions.

Adding the constraints to the L-system

Once defined, the constraints are added to the turtle head Ĥ:

Ĥ = K K K = ε h ρ r Ĥ + ε d ρ r D + n g j=1 ε g, j ρ r Λj + n r i=1 ε r,i ρ i Λi ρ r = n r i=1 (1 -ε r,i ρ i ) (10)
with Ĥ the new turtle head, n g the number of global constraints and n r the number of relative constraints. ε g, j is a user-defined weight for each global constraint j. ε r,i is a userdefined weight for each relative constraint i. ε h and ε d are the user-defined weights for the L-system and the global direction. ρ i is the magnitude of the constraint i and Λi the constraint direction. The magnitude ρ r decreases the turtle's head magnitude, the global direction magnitude and the global constraint magnitudes when a relative constraint has a high magnitude. This improves the conditioning by limiting the influence of the other constraints when getting close to the considered data.

Constraint setting

Constraints call for numerous parameters, which limits their practicality. We have defined default values for the data of section 2.2 [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF]. For most cases, they should be sufficient to ensure the conditioning, and can be adjusted otherwise. Among all the parameters, ten are still required as input (table 3) when constraining all these data. Four parameters are the weights ε of the constraints to adjust their relative importance. The channel global direction is already required by the L-system. The channel minimal and maximal width and minimal thickness come from the distributions required by the L-system. The bandwidth for the well channel data is the only parameter with the constraint weights that is not shared with the L-system. It controls the lateral extension of the perception area and influences the well connections through the channels.

A

The method was implemented in C++ in the Gocad plug-in ConnectO. The NURBS channel model was implemented by Jérémy Ruiu in the Gocad plug-in GoNURBS [START_REF] Ruiu | Modeling channel forms and related sedimentary objects using a boundary representation based on non-uniform rational B-splines[END_REF]. This section develops an application to a synthetic case combining different data. See [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF] for applications to explore parameter and conditioning effects.

Dataset

The synthetic case is build upon submarine channels that migrate within a confining valley. Such channels often display two migration patterns [e.g., [START_REF] Mchargue | Architecture of turbidite channel systems on the continental slope: Patterns and predictions[END_REF]:

• In the lower valley part, channels mostly migrate laterally, with some disorganized stacking between the channels due to discrete migration. This distributes the sandy deposits over the whole valley width.

• In the upper part, vertical migration dominates, with the development of overbanks that confine the channels. The sandy deposits are localized over a smaller section of the valley.

The data set emulates these features through (figure 6):

• A hexahedral grid aligned along its margins, simulated with a L-system and NURBS to materialize the valley.

• A sand probability cube inside the valley grid, with high sand probabilities along the whole valley bottom and more localized at the top. It comes from a sequential Gaussian simulation conditioned to a property whose top derives from a channel simulated with a L-system and whose bottom is evenly filled.

• 234 channel and 391 inter-channel data points extracted along 20 random well paths from a realization conditioned to the sand probability.

The associated constraints are those from table 2. Table A.1 summarizes the input parameters.

Simulation principle

Channels are simulated inside the grid parametric space as long as some well channel data remain unconditioned. Their initial position is drawn close to channel data. The L-system growth is then influenced by the confining valley, the sand probability cube, and the well data, but not by the previously simulated channels. Only an unconditioned channel datum can attract a channel. Once all the channel data are conditioned, channels are still simulated up to a target number of channels -here 40 -if not already reached. The initial position of these last channels is randomly drawn inside the grid. The probability cube influences this draw, so that channels appear preferentially in the high sand probability areas.

Simulation results

The simulated channels remain sinuous despite the significant number of constraints (figure 7). The channel global direction skews the bends, similarly to the observations on real cases (figure 8). However, the bends are deformed depending on the constraint direction, which is not necessarily oriented along the global flow direction. With a high global direction weight, the bends have opposite asymmetries on the two branches (figure 8). For now, no solution has been found to control the bend asymmetry.

The channel centerlines never go outside the valley, but sometimes the channel objects do. This often occurs when an inter-channel datum is too close to a margin: the two constraints may compete with each other, and the channel ends up violating both (figure 9,d). The channels honor all the channel data, as the simulation would not stop otherwise. But the conditioning area is based on a horizontal and a vertical distance: it is rectangular, not channel-shaped (figure 9,c). For the realization in figure 7, 15 over 234 channel data points are outside a channel and would require a post-process to strictly condition them. Similarly, some channels condition 30 over 391 inter-channel data. This is either a failure of the repulsive process (figure 9, b), or a channel that conditions both a channel and an inter-channel datum (figure 9,a). This last case arises because the simulation is meant to favor channel data conditioning: an inter-channel datum below the channel data to condition is ignored.

59,38 channels are simulated on average over 100 realizations, with a standard deviation of 3.64. No realization has less than 49 channels, far from the 40 channels required in input. And the well data come from a realization of 40 channels. The number of channels increases due to well channel data conditioning.

Effect of the constraints

The constraints deform the channel morphology simulated by the L-system to condition the data. E-maps of realizations with different constraints (figure 10) highlight this impact. With the current setting, the margins have a significant repulsive effect on the channels, preventing them from equally covering the entire valley (figure 10,a). This is why the channels tend to gather along the center of the valley bottom, in contradiction with the sand probability cube (figures 6 and 10, b). However, the channels strictly follow the upper trend of the sand probability. The resulting E-map is less blurry than the sand probability, but the channel filling should be simulated for further comparison. Adding well sedimentary data causes the channels to locally divert from the upper trend of the sand probability (figure 10,c). It exemplifies the constraint competition with two data sets that are not fully compatible. Here priority is given to well data, sometimes placing a lower emphasis on the sand probability cube.

A comparison with a conditioning by rejection sampling helps to further analyze the bias introduced by the constraints on the channel morphology (figure 11). The rejection sampling is a simplistic but unbiased conditioning process: an unconditional channel is simulated and kept if it fits the data. Otherwise, a new channel is simulated. Here the channel sections within the inter-well area appear slightly less variable with the constraint than with the rejection sampling. It comes from the straighter morphology of the channels due to the constraint. This bias could be reduced by adapting the constraint setting, especially the magnitude.

Simulation time

Simulating2 channels in the confining valley is quite fast (table 4). Adding constraints impacts the computation time, but it remains acceptable for stochastic purposes. For the probability cube, the main time increase comes from the search for the constraining cell in the perception area. A small area length l a limits the computation time but the conditioning becomes more sensitive to small scale variations. The time increase when adding well data mainly comes from the higher number of simulated channels to condition all the channel data. The NURBS generation is the most time consuming process.

The comparison with the simplistic rejection sampling of fig- ure 11 yet highlights the significant speed of the conditioning process by constraints (table 5). On average, the constraint is 100 times faster than the rejection sampling. And this case has fewer data than the confining valley case, in which the rejection sampling is not an option. The smaller computation times of the NURBS surfaces under the constraint influence come from the straighter channels.

D

The synthetic case illustrates the method ability to simulate channels while conditioning various and numerous data. This section discusses some aspects of the simulation process.

About the use of L-systems

Formal grammars remain uncommon in geological simulations, despite a successful use in other domains. [START_REF] Hill | Describing and generating facies models for reservoir characterisation: 2D map view[END_REF] use a different kind of formal grammar to simulate channels: the plex grammar [START_REF] Feder | Plex languages[END_REF]. In the L-system formalism, a letter has two attaching points, with a letter to its left and one to its right. The plex grammar generalizes this principle to connections between an arbitrary number of letters. It better handles three-dimensional shapes, but the rules quickly become unwieldy. [START_REF] Hill | Describing and generating facies models for reservoir characterisation: 2D map view[END_REF] bypass this issue by using a training model, from which the rules are inferred. Our approach complements it by simulating channels from scratch, avoiding a training model that may be difficult to get. With predefined rules, the L-system itself does not require more parameters than other object-based approaches [e.g., [START_REF] Viseur | Simulation stochastique basée-objet de chenaux[END_REF][START_REF] Deutsch | FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems[END_REF][START_REF] Hassanpour | Improved geostatistical models of inclined heterolithic strata for McMurray Formation, Alberta, Canada[END_REF]. The bend half-wavelength, bend amplitude and channel width can be inferred from seismic data, the channel thickness from well data [e.g., [START_REF] Wonham | 3D sedimentary evolution of a canyon fill (Lower Miocene-age) from the Mandorove Formation, offshore Gabon[END_REF]. When the well data are too sparse or the seismic resolution too low, these parameters can come from seismic or outcrop data of analog settings [e.g., [START_REF] Mchargue | Architecture of turbidite channel systems on the continental slope: Patterns and predictions[END_REF][START_REF] Colombera | A database approach for constraining stochastic simulations of the sedimentary heterogeneity of fluvial reservoirs[END_REF].

One main advantage of L-systems or other formal grammars is the possibility to change the rules. Modifying the channel morphology or introducing non-stationarity along the channel path [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF] becomes easy once the formalism is known. This brings more flexibility to the object definition than in other approaches [e.g., [START_REF] Viseur | Simulation stochastique basée-objet de chenaux[END_REF][START_REF] Deutsch | FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems[END_REF][START_REF] Hassanpour | Improved geostatistical models of inclined heterolithic strata for McMurray Formation, Alberta, Canada[END_REF] and facilitates user interaction. It is one step toward dealing with one major drawback of objectbased methods: the need for predefined geometry for the ob- Asymmetric bends on a seismic section from Angola (modified from [START_REF] Kolla | Evolution of Deep-Water Tertiary Sinuous Channels Offshore Angola (West Africa) and Implications for Reservoir Architecture[END_REF]) and below a channel simulated with a L-system showing opposite bend asymmetries from the initiation point.

jects to simulate. Here focusing on the channel morphology helps to preserve the rule simplicity. But finding simple rules to simulate the whole three-dimensional shape of the channel could be a step forward to a more interactive modeling process, such as what is already done in biology [e.g., [START_REF] Prusinkiewicz | The Algorithmic Beauty of Plants[END_REF][START_REF] Prusinkiewicz | The use of positional information in the modeling of plants[END_REF][START_REF] Longay | Treesketch: interactive procedural modeling of trees on a tablet[END_REF]. Although channels simulated with a L-system can be used to initiate a migration process [START_REF] Rongier | A geostatistical approach to the simulation of stacked channels[END_REF], it should be possible to define rules that simulate an entire sequence of migrating channels instead of a single channel.

In both cases, the resulting sedimentary deposits could be deduced from channel intersections [START_REF] Ruiu | Modeling channel forms and related sedimentary objects using a boundary representation based on non-uniform rational B-splines[END_REF].

About constraints and conditioning

Conditioning through constraints relies on a sum of vectors. Thus, it is a simple and fast process. Getting the constraining vector is the only aspect that may slow down the process. The concept of constraining element and vector is flexible enough to handle a wide variety of data. Here we introduce the conditioning of well channel and inter-channel data, probability cubes and confining structures. Other elements can be taken into account, such as well static connectivity data [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF]. Now it should be extended to the conditioning of partially interpreted channels from seismic data or other well sedimentary data, such as levee and lobe. Adding all these data can help to reproduce the channel architecture, as shown in the case study, and get closer to the channel stacking deriving from channel migration. To better reproduce this stacking, the previously simulated channels could become attractive constraints for the newly simulated channels. Constraints could then reproduce the result of channel migration and avulsion without modifying the L-system rules.

But this flexibility calls for numerous parameters, with some not necessarily easy to determine, such as the perception area. The solution proposed here is to predefine the constraints and their parameters. But depending on the studied case and the number of constraints, these parameters may not be compatible enough, leading to a poor conditioning, and should be modified. All this requires more work to reduce the number of constraint parameters or to infer some of them automatically.

A similar conditioning process is already applied to channel simulation [START_REF] Lopez | Modélisation de réservoirs chenalisés méandriformes : une approche génétique et stochastique[END_REF][START_REF] Pyrcz | ALLUVSIM: A program for eventbased stochastic modeling of fluvial depositional systems[END_REF] through a lateral deformation of the final channel. Here the constraints are directly oriented toward the data as the channel grows: the channel can go straightly from one data to the other to ensure the conditioning. Compared to rejection sampling [e.g., Deutsch andTran, 2002, Hassanpour et al., 2013], constraints are faster and more flexible, especially if the specified channel morphology is inconsistent with the data. More recent sampling approaches manage to condition high density well data [START_REF] Hauge | Object-Based Modeling with Dense Well Data[END_REF]. Although using object-based instead of cell-based methods is questionable with a densely-drilled fieldwhen cell-based methods correctly reproduce the channel continuity thanks to the data density -constraints have also shown promising results with denser well data than in our case study [START_REF] Rongier | Connectivity of channelized sedimentary bodies: analysis and simulation strategies in subsurface modeling[END_REF]. And constraints condition data other than wells more easily compared to other object-based approaches. Compared to multiple-point simulation [e.g., [START_REF] Strebelle | Conditional simulation of complex geological structures using multiple-point statistics[END_REF][START_REF] Arpat | A Multiple-scale, Pattern-based Approach to Sequential Simulation[END_REF][START_REF] Mariethoz | The Direct Sampling method to perform multiple-point geostatistical simulations[END_REF], constraints ensure the preservation of the channel continuity and shape.

Improving the conditioning

The L-system is not always able to condition all the well data, and a channel may condition an inter-channel data or miss a channel data (figure 9.b and c). If some improvements are undoubtedly possible, other solutions exist to ensure data conditioning. A simple solution is to re-simulate the channel. This is worth considering thanks to the low simulation time, but only for significant conditioning issues. Another possibility is to re-simulate the channel width or thickness. This is useful if a channel goes slightly outside its confinement for instance. Finally, the last possibility is to rely on the NURBS, which are very easy to deform. Such deformation process could be useful to condition other data, such as channel orientations along a well.

The sand probability cube derived from seismic data is not equivalent to a channel probability cube. Indeed, sedimentary bodies are heterogeneously filled: levees may contain sandy deposits and channels muddy deposits. Channel fill and sedimentary bodies around channels should also be simulated to further study the sand probability reproduction. It implies that the channel E-map may only reproduce the global trend of the probability cube. If the probability cube is composed of discontinuous patches of high probabilities, simulating channel objects may not be the best option. Other solutions should be considered, for instance simulating lateral accretion packages [START_REF] Hill | Describing and generating facies models for reservoir characterisation: 2D map view[END_REF]Griffiths, 2009, Hassanpour et al., 2013]. In any case, making sure that the local variations of that probability are reproduced remains here an open issue.

Another improvement concerns the stopping criterion. For now, conditioning well channel data leads to more channels than required. Improving the conditioning can help to stay closer to the stopping criterion. Conditioning a global channel proportion or a net-to-gross would be more appropriate due to their impact on flow simulation. For this, other object-based methods rely on iterative processes [e.g., Viseur, 2001, Deutsch and[START_REF] Deutsch | FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems[END_REF].

C

This work introduces a new method for channel simulation through a formal grammar. The formal grammar, called the Lindenmayer system, simulates the channel morphology using predefined rules. Modifying the rules transforms the morphology, which brings new flexibility to object simulation. Attractive and repulsive constraints ensure the conditioning by influencing the channel development. They call for numerous parameters, but this approach enables to condition various data within the same framework. Most of the parameters can be predefined for an easier use. The conditioning deforms the channel morphology and a slight statistical bias. But it ensures the conditioning: if the priority is given to the data constraint, the channel straightly goes to the data or avoids it depending on the constraint type. The channel sinuosity can be lost, but it is less essential in flow simulation than the channel continuity.

Many improvements are still possible. Some work should be done around NURBS deformation to perfectly fit the data at small scale. And the conditioning process could be improved to less deform the channel morphology. From this point of view, a lot can be borrowed from robot motion planning, as the principles are pretty similar. It also includes introducing other data, such as the net-to-gross. Other sedimentary structures should be integrated within the L-system rules, such as lobes or levees, whose NURBS parameterizations already exist [START_REF] Ruiu | Modeling channel forms and related sedimentary objects using a boundary representation based on non-uniform rational B-splines[END_REF].

A

This work was performed in the frame of the RING project at Université de Lorraine. We would like to thank the industrial and academic sponsors of the Gocad Research Consortium managed by ASGA for their support and Paradigm for providing the SKUA-GOCAD software and API. We would also like to thank the reviewers for their constructive comments which helped improve this paper.

A A S

Table A.1 Parameters to simulate the channels of the synthetic case. is a triangular distribution with a minimum, a mode and a maximum.

Simulation parameters Values

Number of channels 40

L-system

Global direction (in °) 90

Global direction weight 0.25 Default segment length (in cell) 6

Half-wavelength (in cell) (10,15,25) Amplitude (in cell) (0,4,7)

Deviation angle (in °) (0,0.57,5.7)

L-system weight 1

Channel sections

Channel width (in cell) (5,6,8) Channel width range (in cell) (10,15,20) Curvature weight 0.75

Channel thickness (in cell) (1.5,2,2.5) Channel thickness range (in cell)

Curvature weight 0.75

Asymmetry aspect ratio 0.5

Constraints

Confinement weight 1

Channel data weight 1

Channel data bandwidth (in cell) 30

Inter-channel data weight 1

Sand probability weight 1

Figure 1

 1 Figure 1Channel discretization for L-system simulation. The red dots locate the channel sections obtained from the interpretation. l s is the distance between two successive sections, i.e., the channel segment length. α is the angle between two successive sections, i.e., the change of orientation between two successive channel segments. The other parameters characterize a channel bend, with l B the bend length, λ the bend half-wavelength, ∆ the bend amplitude, and r = 1/c the radius of curvature -radius of a circle fitted at the bend apex -with c the curvature.

  Figure 3Constraint addition to a L-system, with Ĥ the current turtle's head, Ĥ the new turtle's head influenced by the constraint, D the channel global direction, Λatt the direction for an attractive constraint, deduced from V , Λrep the direction for a repulsive constraint, deduced from -V and D, and V the vector between the turtle's head and the constraining element.

Figure 4 Figure 5

 45 Figure 4Parameters defining the perception area, with d its direction, σ its direction tolerance, w a its maximal half-width, l a its length and t a its thickness. Ĥ is the turtle's head.Inner areaTransition area Outer area

Figure 6

 6 Figure 6 Data set of the application: a curvilinear grid representing a confining valley with a sand probability cube and twenty wells.

Figure 7

 7 Figure 7Realization conditioned to the sand probability cube and well sedimentary data within the confining valley The highest channel -on the right -must remain inside the valley while following the high sand probability without conditioning the 14 inter-channel data and with 6 channel data that it may condition. The two channel data outside this channel condition another channel.

Figure 8

 8 Figure 8Asymmetric bends on a seismic section from Angola (modified from[START_REF] Kolla | Evolution of Deep-Water Tertiary Sinuous Channels Offshore Angola (West Africa) and Implications for Reservoir Architecture[END_REF]) and below a channel simulated with a L-system showing opposite bend asymmetries from the initiation point.

Figure 9 Figure 10 Figure 11 E

 91011 Figure 9Details on well conditioning. All the channels are not showed for clarity. a. The four channels at the top perfectly condition the channel data. In the middle, the last channel condition both a channel data and an inter-channel data: currently the process can not always prevent the inter-channel data conditioning. b. Five channels are close to a set of inter-channel data. One of them condition the inter-channel data, showing a failure in the repulsion process. c. The channel conditions the two channel data as defined in the method. However, the channel data at the bottom is not within the channel and is not stricto sensu conditioned. d. Example of incompatibility between the L-system path and the data.

Table 1 L

 1 -system alphabet for channel simulation.

Table 2

 2 Constraints used in section 3.

	Constraint	Type	Element
	Well channel data	Relative & attractive Well point
	Well inter-channel data	Relative & repulsive Well point
	Probability cube	Global & attractive	Grid cell
	Confinement (margin 1) Relative & repulsive Grid border
	Confinement (margin 2) Relative & repulsive Grid border

Table 3

 3 Input parameters for the constraints of section 2.2, with ε i the constraints weights, t min the minimal channel thickness, D the channel global direction, w a,C D the perception area bandwidth for the well channel data and w min and w max the minimal and maximal channel widths.

			Input parameters
	Well channel data	ε C D	D w a,C D w min t min
	Well inter-channel data	ε I C D D -	w max t min
	Probability cube	ε P C	D -	-	t min
	Confinement (margin 1) ε C	D -	w max -
	Confinement (margin 2) ε C	D -	w max -

Table 4

 4 Simulation time for a realization in the confining valley with different constraints. The given values are the average time over 100 realizations and the standard deviation.

	Constraints	L-system simulation (in s) NURBS generation (in s)
	Confinement only	3.17 ± 0.15	16.24 ± 0.22
	Confinement and probability cube	4.73 ± 0.10	14.57 ± 0.15
	Confinement, probability cube and well data 6.64 ± 0.45	21.47 ± 1.31

Table 5

 5 Simulation time for a realization containing one channel with two different conditioning processes. The given values are the average time over 1 000 realizations and the standard deviation.

Full text available at: https://hal.archives-ouvertes.fr/tel-01371350

On a 64-bit Linux system with a 2.10 GHz processor Intel® Core™ i7-3612QM and 6 GB of RAM.


Jean Ogier

Du Terrail

Frederic Jurie
email: frederic.jurie@unicaen.fr

ON THE USE OF DEEP NEURAL NETWORKS FOR THE DETECTION OF SMALL VEHICLES IN ORTHO-IMAGES

Keywords: target detection, deep-learning, aerialimagery

. This work supports an extensive study of the problems one might face when applying deep neural networks with low resolution and scarce data and proposes some solutions. One of the contributions of this paper is a network severely outperforming the state-of-the-art while being much simpler to implement and a lot faster than competitive approaches. We also list the limitations of this approach and provide several new ideas to further improve our results.

INTRODUCTION AND RELATED WORK

Although automatic detection of objects in images is an old problem and benchmarks have been around for a while [START_REF] Everingham | The PASCAL Visual Object Classes Challenge 2012[END_REF], the computer vision community only recently began to put the focus on detecting small objects. Modern benchmarks like [START_REF] Lin | Microsoft COCO: Common Objects in Context[END_REF] that favored small objects have proven to be surprisingly challenging and to this day no approach has clearly gained the upper-hand on the others. In spite of the progress of the technologies behind satellite optics it is likely that aerial imagery pictures will always be filled with small objects with relatively low resolution. Furthermore, as aforementioned, the vehicles to detect occupy a few hundreds pixels on often very large images (around a million pixels) so most pixels/windows on such images will have to be classified as belonging to a sink class that we will call background. Those very peculiar conditions explain why one cannot directly apply winning methods on [START_REF] Everingham | The PASCAL Visual Object Classes Challenge 2012[END_REF] or [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF] like [START_REF] Shaoqing Ren | Faster R-CNN: towards real-time object detection with region proposal networks[END_REF] on such benchmarks.

Modern pipelines for detection can be roughly classified in three main groups. Grid-based regression methods: the user defines a grid on the image and regress bounding boxes based on offsets w.r.t this grid [START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF][START_REF] Szegedy | Scalable, High-Quality Object Detection[END_REF][START_REF] Liu | SSD: single shot multibox detector[END_REF]. Region proposals based methods: most of them consisting of a cascade starting with a class agnostic region proposals step [START_REF] Cheng | BING: Binarized normed gradients for objectness estimation at 300fps[END_REF][START_REF] Uijlings | Selective search for object recognition[END_REF][START_REF] Zitnick | Edge boxes: Locating object proposals from edges[END_REF] followed by a classifier [START_REF] Girshick | Rich feature hierarchies for accurate object detection and semantic segmentation[END_REF][START_REF] He | Spatial pyramid pooling in deep convolutional networks for visual recognition[END_REF][START_REF] Girshick | Fast R-CNN[END_REF]. This category also includes the recent development that began with [START_REF] Shaoqing Ren | Faster R-CNN: towards real-time object detection with region proposal networks[END_REF] and that consists of learning the region proposal part together with the classification (see e.g. [START_REF] Dai | R-FCN: object detection via region-based fully convolutional networks[END_REF]). Sliding window based methods [START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF][START_REF] Szegedy | Deep neural networks for object detection[END_REF]. We argue that none of the methods above could be applied without significant modifications to the task of detecting small targets in large images. For the first class of methods the grid used would have to be excessively large with a closeknit network, which would require us to work on pooling or up sampling the final output and that goes against the philosophy of grid-based methods, which is to be coarse to be fast. In the second class of methods, problem specific region proposal algorithms would have to be designed so that they could detect small low-resolution objects with high recall. For instance, we think that the work we propose could be used as a first step in such pipelines.

Vehicle detection on aerial imagery has been extensively studied in [START_REF] Zhao | Car detection in low resolution aerial images[END_REF][START_REF] Leberl | Recognizing cars in aerial imagery to improve orthophotos[END_REF][START_REF] Eikvil | Classificationbased vehicle detection in high-resolution satellite images[END_REF][START_REF] Kembhavi | Vehicle detection using partial least squares[END_REF][START_REF] Ali | A real-time deformable detector[END_REF]. The most recent articles on the subject are [START_REF] Razakarivony | Apprentissage de variétés pour la Détection et Reconnaissance de véhicules faiblement résolus en imagerie aérienne[END_REF][START_REF] Gleason | Vehicle detection from aerial imagery[END_REF], based on handcrafted descriptors, and [START_REF] Chen | Vehicle detection in satellite images by parallel deep convolutional neural networks[END_REF][START_REF] Mundhenk | A large contextual dataset for classification, detection and counting of cars with deep learning[END_REF][START_REF] Audebert | On the usability of deep networks for object-based image analysis[END_REF] using convolutional neural networks (CNNs). [START_REF] Chen | Vehicle detection in satellite images by parallel deep convolutional neural networks[END_REF] also uses a CNN with a similar architecture, but does not provide any strategy to deal with class imbalance. It uses a different and outdated cost function and is much slower. [START_REF] Mundhenk | A large contextual dataset for classification, detection and counting of cars with deep learning[END_REF] introduces a new and large dataset and reports good performances but the context-based method it used is impractical in our case, furthermore, the metric used to measure perfor-mances is much less demanding in terms of precision of the localization.

This work focusses on developing hard-mining strategies to deal with classes with few examples and is a lot faster than all former approaches due to the small sized network and the fully convolutional inference and is specifically tailored for the detection of small objects.

METHODOLOGY

Introduction

In contrast with all the previous detection results on VEDAI, we do not use any sliding a window over the image to do inference but get rid of this expensive step using fully convolutional networks (FCN), as proposed by [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF]. In terms of architectures, we have experimented with several variants around the simple LeNet-5 architecture [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF]. This architecture seems indeed a good candidate for our detection task as the LeNet5 network has only 2 max-pooling layers (nonoverlapping). Consequently, when used as a FCN the resolution of the output heat maps is 4 times smaller than the original image. In practice there is no need to up sample the heat map by unpooling with max-pooling switches or using transposed convolutions like in [START_REF] Shelhamer | Fully convolutional networks for semantic segmentation[END_REF] or [START_REF] Badrinarayanan | Segnetme: A deep convolutional encoder-decoder architecture for image segmentation[END_REF], nor to use dilated convolutions as in [START_REF] Yu | Multi-scale context aggregation by dilated convolutions[END_REF] nor even to use the shift-and-stitch trick of [START_REF] Sermanet | Overfeat: Integrated recognition, localization and detection using convolutional networks[END_REF]. One key difference of detection on ortho-images and detection on commonly used datasets like MS COCO, is that, because the distance between the sensor/camera and the ground is known in advance, all targets of a same class share the same size in pixels, approximately, and this size can be estimated accurately. Therefore there is no need to adopt a multi-scale approach nor to regress the width and height of the bounding boxes of the vehicles. For this reason, we can set once and for all the size of the extracted patches. By doing that we also set the size of the first fully connected layer. The imbalance between the sink class and the vehicle class being so overwhelming we had to use a multi-step approach to get rid of most of the backgrounds and simplify the classification as noticed by [START_REF] Audebert | On the usability of deep networks for object-based image analysis[END_REF] once the vehicles and the background are identified classification between the different classes is more straightforward.

Training

We first extract positive patches around the targets. The negative patches used are sampled uniformly from the images until we reach a ratio of 5 to 1 (which has been determined by cross-validation). Too many negative samples and the imbalance is too strong, too little and the heterogeneity of the background class is not fully captured. A critical point that is also investigated in section 5 is the imbalance factor applied to the cross-entropy cost. Then the network is trained using stochastic gradient descent (SGD), with a dropout of 0.5 [START_REF] Hinton | Improving neural networks by preventing co-adaptation of feature detectors[END_REF] in the fully connected layers. When it reaches convergence, after around 100 epochs, then we do inference on the images on the training set in a fully convolutional fashion. From 1064x1064 pixels image (images have been up-sampled to 1024x1024 and padded) we get 256x256 heat maps. To select the maximums we do traditional Non-Maximum Suppression (NMS). We then consider the remaining maximum values to be the result of our detection and we evaluate them on the ground truth.

In the process we sample negative patches that have been misclassified and weak true positives. In fact we experimented with 4 different strategies of hard-mining results are presented in section 5. We repeat the process multiple times.

Hard-Mining strategies

Bootstrapping offers a lot of liberties on how the hard examples are chosen. One could for instance pick a limited number of false positives per image or one could fix a threshold and only pick a false positive if its score is superior to a fixed threshold (0.5 for instance). There is also the question of weak true positives whether to include them at all and if yes from which threshold should we pick them. The previous experiments used the thresholds (0.5,0.5) while limiting the number of examples in each image to be 25 at max. We tested the 4 following different hard-mining variants: Strategy 1: the number of hard examples per image is set to be exactly N (with N=25 in our experiments) whatever their scores. No weak true positive are added. Strategy 2: Same strategy that 1 but with weak true positive examples whose scores are less than 0.5 (scores are normalized probabilities). Strategy 3: We chose only hard negative with scores superior to 0.5 and weak true positives with scores less than 0.5. Strategy 4: Same as strategy 3 without weak true positives. These 4 strategies are evaluated in Section 3.

Addressing class imbalance issues

One of the issues with detection is the imbalance between the few present vehicles (targets) on images and the large variety of backgrounds. This issue is even more present in aerial imagery because of the relative size of full images and small vehicles like cars. In this context, VeDAI is especially challenging as each image only contains a few vehicles. There are many techniques to effectively fight this. Controlling vehicles/background proportion in each batch, minority oversampling like [START_REF] Nitesh | Smote: synthetic minority oversampling technique[END_REF], or modifying the cost itself to give more weights to the misclassification of classes that are less present. We chose the latter because of its simplicity. The cost function implemented is the weighted cross-entropy:

L(w) = N i=1 C * t i * log(f w (x i )) + (1 -t i ) * log(f w (x i ))
where x i is one of the N image patches in the training set, f w (x) is the score given by the convolutional neural network to a patch x, t i is 1 when x i belongs to the class under test 0 otherwise, C is a scalar defined as the ratio of Negative (background) vs Positive (vehicles) examples in the batch.

Classification study: rotation invariance, etc.

This part aims at finding the weaknesses of our approach and further motivates a follow-up work on using this baseline network as the first part of a cascade. We created numerous classification sets using collected hard-negative samples. We looked at every possible combination of: contrast normalizing the patches (or not), angle normalizing the targets (meaning all targets have the same orientation)(or not) and shifting the patches from the targets (4 pixels apart in checkerboard distance)(or not). For simplicity we adopted the following notation: S and ¡ S mean respectively with shifted target normalization or without, A and A means angle-normalized patches or all rotations included, and, C and C contrast-normalized patches or no contrast normalization used. These different options are experimented in Section 3

EXPERIMENTS

The VeDAI dataset The VeDAI dataset [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] consists of 1200 images that come in two different resolutions 512x512 and 1024x1024. Every image is available in two versions either colored (RGB) or infrared. All experiments in the paper were conducted on the infrared version of the 512x512 images. The dataset is provided with 10 folds and 1340 cars in total, with an average of around 140 cars to detect per fold. The definition of what a positive detection is pretty different from the standard Intersection Over Union (IOU) criterion adopted by the Pascal VOC or MS COCO. A detection is considered correct if it lies within the ellipse centered on the ground truth and lying inside the edges of the target car (multiple hits on the same targets are counted as False Positives). We are interested in mainly two metrics namely the mAP which is measured across all folds and the recall at low FPPI (recall for a given rate of False Positive Per Image).

Results on VeDAI

We first report general results on VeDAI using our architecture inspired by the LeNet-5 network. On overall, the network contains 2 convolutional layers followed by three fully connected like layers, as detailed Table 1.

All convolutions are VALID type convolutions (no padding) and all fully connected layers are implemented as 1x1 convolutions. In order to make the results as good as possible, we experimented with the following parameters : (i) size of the receptive field, (ii) amount of regularization (iii) depth of each layer. We selected optimal parameters on a validation set. We split the 10 folds provided with the dataset into ten subsets of equal length and chose our parameters by training the network on 9 of them and validating on the tenth.

The selected network reaches a mAP of 77.8±3.3 (see Fig. 1) which is a very great improvement (+12%) from the 2 lists all the published results on VeDAI and ours. We can see that in terms of Average Precision or recall the advantages of our method. [START_REF] Audebert | On the usability of deep networks for object-based image analysis[END_REF] do not report any detection results only classification so we cannot compare to it directly.

Evaluation of the 4 hard-negative-mining strategies As said before, bootstrapping offers a lot of liberties on how the hard examples are chosen. We proposed in Section 2 4 different hard-mining strategies we are now going to evaluate. Table 3 gives the Recall at 0.01 False Positive per Image (FPPI) for different rounds of hard mining. From this table we can make three observations: first, the best overall performance is obtained with S2 and should be preferred. Second, if one wants to limit the number of rounds of hard example mining, he should prefer S3 which gives 0.35 after only 2 rounds. Third, without hard mining, the performance of the network method AP Recall@0.01FPPI DPM [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] 60.5 ± 4.2 13.4 ± 6.8 SVM+HOG31 [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] 55.4 ± 2.6 7.8 ± 5.5 SVM+LBP [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] 51.7 ± 5.2 5.5 ± 2.2 SVM+LTP [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] 60.4 ± 4.0 9.3 ± 3.7 SVM+HOG31+LBP [START_REF] Razakarivony | Vehicle Detection in Aerial Imagery: A small target detection benchmark[END_REF] 61.3 ± 3.9 8.3 ± 5.2 SVM Fusion AED (HOG) [START_REF] Razakarivony | Apprentissage de varietes pour la Detection et Reconnaissance de vehicules faiblement resolus en imagerie aerienne[END_REF] Compensating the translations, rotations, etc.. We experimented the 3 normalization strategies given Section 2.5. As said before, to make the evaluation of the performance faster, we experimented these alternatives on a classification task. For building this classification set, we used all the targets from fold 1 and added hard negative examples. There is 50000 examples in the training set with a ratio of negative over positive of 70 it is 5 times more than in the course of hard-mining passes as we had to sample only one positive example per target present instead of 5 as we have when we start detection training. We measure the classic classification accuracy but as the number of backgrounds is 70 times the number of targets, it is actually not very informative. It is more helpful to have also the accuracy on the target class patches, which verifies that the classifier would not get away with classifying everything as background and still get a good accuracy. We added the recall for 0.001 FP meaning the recall of the target classes patch obtained for a given proportion of 0.001 false positives in the validation set for the same reason. The results are given in Table 4. It is also interesting to see what would the network that was trained on perfectly centered patches do on shifted targets (it often happens when the equivalent sliding window step is not 1). We get 16.

CONCLUSIONS

We have presented state-of-the-art results on a challenging benchmark and insights on how to tackle detection in the difficult context of few examples learning and extreme imbalance. This study cuts across the many articles using transfer learning to avoid having to deal with training a network from scratch. We are also confident that having found in Section3 that angle normalization has such a strong effect on performances we could improve a lot on these results by using a cascade.

Fig. 1 .

 1 Fig. 1. Precision-Recall curves given by our RPN network on the 10 folds of the VEDAI dataset. The mAP is of 0.78± 0.03.

Table 1 .

 1 Architecture of our network, inspired by LeNet-5

	Name	Type	Filter Size	Input	Output
			Stride	Size	Output
	Conv1	convolution	5x5/1	45x45x3	41x41x96
	Pool1	max-pooling	2x2/2	41x41x96	20x20x96
	Conv2	convolution	5x5/1	20x20x96	16x16x192
	Pool2	max-pooling	2x2/2	16x16x192	8x8x192
	fc1	fully connected	8x8/1	8x8x192	1x1x384
	fc2	fully connected	1x1/1	1x1x384	1x1x84
	fc3	fully connected	1x1/1	1x1x84	1x1x2
	Softmax	Softmax	None	1x1x2	1x1x2

Table 2 .

 2 Comparisons with related works

					69.6 ± 3.4	20.4 ± 6.2
		Ours			77.80 ± 3.3	31.04 ± 11
	Strategy R1	R2	R3	R4	R5	R6
	S1	0.00 0.21 0.39 0.36 0.28 0.25
	S2	0.00 0.10 0.17 0.43 0.27 0.32
	S3	0.00 0.35 0.34 0.31 0.25 0.29
	S4	0.00 0.23 0.32 0.23 0.32 0.32

Table 3 .

 3 Recall@0.01FPPI w.r.t. the number of passes for the 4 strategies, on a validation set is very low.

Table 4 .

 4 Learning variant/invariant representations. See Section 3 for details. are centered and normalized in rotations. We ran additional experiments to see if the angle normalization step of the cascade was really necessary so we added a data-augmentation module which rotated the patches in the current batch by a random angle uniformly chosen between 0 and 360. The network that performed poorly (53.91 accuracy on positive examples and 28.02 recall at 0.001 FPPI) gave reasonable results with 88.77 accuracy on the positive examples. So if we can center and normalize detected results we could use this trained classifier on top and it would better the results.

	5 accuracy on the positive examples
	(instead of 53.91) and a recall at 0.001 FPPI of 18.20 (instead
	of 28.02). There is a severe degradation of the performances.
	The same phenomenon is observed when we train on all ori-
	entations and test on angle normalized patches. Accuracy on
	positive examples drops to 52.67 (instead of 82.71). The re-
	call at 0.001 FPPI is also impacted 55.8 (instead of 63.45).
	Therefore this network performs much better when the targets

* This work was supported by ANRT through the CIFRE sponsorship No 2015/1008 and by SAFRAN Electronics and Defense.