Szymon Racewicz
email: szymon.racewicz@g2elab.inpg.fr

Delphine Riu
email: delphine.riu@g2elab.inpg.fr

Nicolas Retiere

Gérard Meunier

Patrice Labie

Piotr J Chrzan

HALF-ORDER ANALYTICAL MODELLING OF FERROMAGNETIC SHEET. VALIDATION WITH FINITE ELEMENT SIMULATION

This paper presents a half-order modelling approach applied to a saturated ferromagnetic sheet. Several analytical models are then developed and compared for 1D formulation. In order to identify parameters and valid such models, the saturated ferromagnetic sheet has been simulated using Flux® finite element simulation software.

I. INTRODUCTION

Nonlinearity is a phenomenon which has to be taken into account in diverse applications, e.g. design of electric machines or transformers [START_REF] Mayergoyz | Nonlinear diffusion of electromagnetic fields: with applications to eddy currents and superconductivity[END_REF]. Then, this phenomenon cannot be neglected in order to obtain high precision models. Non integer order systems are used to model diffusion in conductive parts of electrical machines as they lead to more compact and knowledge models but also to improve their precision [START_REF] Riu | Induced currents modeling by half-order systems application to hydro-and turbo-alternators[END_REF]. Nevertheless non integer order models rest valid for linear cases. It is then interesting to study the influence saturation on such model parameters.

In this article, we first present the linear half-order model of a ferromagnetic sheet deduced from the diffusion equation of the magnetic field. Then, a new formulation taken into account the skin effect and the saturation of the magnetic materials is presented and compared to classical analytical models. Finally, the non linear halforder sheet model is identified and validated from a finite elements simulation.

The studied ferromagnetic sheet is presented in Figure 1. It is characterized by permeability (µ), resistivity (ρ) and surrounded by a n-turn winding.

     ⋅ ⋅ ⋅ = ⋅ ⋅ ⋅ = e j L l e n L µ ρ ω α µ r (2)
L 0 is the DC inductance of the sheet; an integer order approximation of this impedance from serial decomposition of tanh function can then be used in simulation software.

II.2. Half-order model

From a second-order limited development of tanh function, another approximation of the analytical impedance is presented introducing half-order systems [START_REF] Racewicz | Half-order modelling of turboalternators -An adapted method of parameter identification[END_REF].

2 0 0 0 2 / 1 2 1 1 with 1       ⋅ ⋅ = ⋅ + ⋅ = e j j L Z µ ρ ω ω ω ω r (3)
Figure 2 presents a comparison between the analytical expression and its approximation using half-order system. The half-order model is then a good approximation of the analytical result along a wide frequency range, even if a small gap around the cut-off frequency can be noticed.

. State of art on 1D non-linear models

There are several ways to model impedance in saturated conditions. First of all, we consider that the first harmonic hypothesis is available; therefore, thanks to the method of equivalent permeability used by simulation software, all electrical variables are considered to be sinusoidal.

Classical impedance variations with current and frequency can be by equation ( 4) :

ω ω σ j l i k L j i L i Z n static ⋅       + ⋅ + = ⋅ = 1 ' ) ( ) , ( 0 r (4)
where: L 0 ', l σ , k and n are parameters identified from measurements [START_REF] Keyhani | Igspice simulations of induction machines with saturable inductances[END_REF]. It is also possible to obtain non-linear impedance model by resolving with severe approximations the non-linear diffusion equation of the electromagnetic field [START_REF] Mayergoyz | Nonlinear diffusion of electromagnetic fields: with applications to eddy currents and superconductivity[END_REF]. Nevertheless, these analytical models are too simple because of the approximative saturation curve. They only allow to characterise the decrease of the DC inductance with current increase.

III.2. Non-linear half-order model

The authors propose a new analytical model of the impedance taking into account frequency phenomena in non-linear conditions. Then, the structure of the presented in equation (3) remains identical, but parameters L 0 and ω 0 vary with the field current. This model is then described by:

) ( 1 ) ( ) , ( 0 0 2 / 1 i j j i L i Z ω ω ω ω ⋅ + ⋅ = r (5)
We suppose that: L 0 (i)=L static (i) with l σ =0 and L 0 '=L 0 . Moreover L 0 and ω 0 are supposed to be linked as ω 0 (i)=m/L 0 (i), where m is a parameter to identify (Fig. 5). Finally, the half-order model of saturated impedance is given by (6). It is characterized by three parameters to identify: k, n and m.

) 1 ( ' 1 1 ' ) , ( 0 0 2 / 1 n n i k m L j j i k L i Z ⋅ + ⋅ ⋅ + ⋅ ⋅ + = ω ω ω r (6)

III.3. Validation with finite element calculations

In order to validate the new model of impedance, finite element simulations were done (Fig. 3). The impedances were calculated over the range of frequency from 0.01 Hz 1000 Hz. The current flowing in a one turn winding varies from 0.01 A to 20 A, which corresponds to a magnetic field between 0.05 and 100 A/m. In Figure 4, the simulation results are shown. Using these data, the set of parameters (k, n, m) of the new model presented in (6) have been identified for one chosen frequency equal to 50 Hz. The identified parameters are given in Table 1. Comparison of the identified model for whole range of frequency with finite element simulations is shown in Figure 6 for two values of field current (0.01 and 20 A). It can be noticed that the model seems to agree with simulations. Then, the asymptote of the impedance is still equal to 10 dB per decade. This model will be soon included in the equivalent circuits of a synchronous machine taking into account on the one hand eddycurrents in massive parts of the rotor or rotor bars and, on the other hand, saturation of the magnetic field.

Fig. 1 .

 1 Fig.1. Studied ferromagnetic sheet II. IMPEDANCE MODELLING IN LINEAR CONDITIONS II.1. Classical analytical models Classical analytical model of impedance can be expressed by equation (1) [3].

  ] P h a s e [d e g ] Classical model Half-order model

Fig. 2 .

 2 Fig.2. Classical and half-order models of the sheet III. NON-LINEAR AND NON-INTEGER MODELLING III.1. State of art on 1D non-linear models

Fig. 3 .

 3 Fig.3. Modelled ferromagnetic sheet with its saturation curve

Fig. 4 .Fig. 5 .Fig. 6 .

 456 Fig.4. First-harmonic model of the saturated impedance


Benjamin Apffel

Filip Novkoski

Emmanuel Fort
email: emmanuel.fort@espci.fr

Floating under a levitating liquid

When placed upside down, a liquid surface is known to fall if it is above a certain size. Gravity acting on the lower liquid interface has a destabilizing effect called Rayleigh-Taylor instability 1,2 .

 or stabilization of the liquid layer above sinking objects 20 . Here, we show that the dynamic stabilization creates a stable buoyancy on the lower interface as if gravity were inverted. Bodies can thus float upside down on

the lower interface of levitating liquid layers. We take advantage of the excitation resonance of the supporting air layer to perform experiments with large levitating liquid volumes up to half a liter and a width up to 20 cm. The experimental findings are explained by the dynamical stabilization of the buoyant equilibrium which exists on the lower interface. This model also predicts a minimal excitation needed to withstand falling which is dependent on the floater mass. Experimental observations confirm the possibility of selective falling of heavy bodies. Our findings invite us to rethink all interfacial phenomena in this exotic and counter-intuitive stable configuration.

Maintaining a liquid upside-down is challenging but various situations in which the liquid can be sustained are known. In the case of a limited surface size, capillary forces have a stabilizing effect [START_REF] Young | An essay on the cohesion of fluids[END_REF][START_REF] Thomson | On the equilibrium of vapour at a curved surface of liquid[END_REF][START_REF] Gennes | Capillarity and Wetting Phenomena: Drops, Bubbles, Pearls, Waves[END_REF] . Alternatively, for large but thin liquid layers suspended under a plate, capillarity opposes gravity [START_REF] Fermigier | Two-dimensional patterns in Rayleigh-Taylor instability of a thin layer[END_REF][START_REF] Myshkis | Mathematical Theory of Capillary Phenomena[END_REF] . In the latter case, the liquid interface does not stay flat but is destabilized in a regular pattern of hanging droplets. This instability driven by gravity, known as the Rayleigh-Taylor instability, occurs at the interface between two fluids whenever a denser one is placed over a lighter one [START_REF] Rayleigh | Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density[END_REF][START_REF] Lewis | The instability of liquid surfaces when accelerated in a direction perpendicular to their planes[END_REF] . Several approaches have been used to stabilize the liquid layer such as temperature gradients [START_REF] Burgess | Suppression of Dripping from a Ceiling[END_REF] , electric [START_REF] Petropoulos | On the control and suppression of the Rayleigh-Taylor instability using electric fields[END_REF] or magnetic fields [START_REF] Rannacher | Suppressing the Rayleigh-Taylor instability with a rotating magnetic field[END_REF] , rotational motion [START_REF] Tao | Nonlinear Rayleigh-Taylor instability of rotating inviscid fluids[END_REF] and vertical vibrations [START_REF] Wolf | The dynamic stabilization of the Rayleigh-Taylor instability and the corresponding dynamic equilibrium[END_REF][START_REF] Wolf | Dynamic Stabilization of the Interchange Instability of a Liquid-Gas Interface[END_REF][START_REF] Lapuerta | Control of Rayleigh-Taylor instability by vertical vibration in large aspect ratio containers[END_REF][START_REF] Kumar | Mechanism for the Faraday instability in viscous liquids[END_REF][START_REF] Pototsky | Faraday instability of a two-layer liquid film with a free upper surface[END_REF][START_REF] Pototsky | Vibration-induced floatation of a heavy liquid drop on a lighter liquid film[END_REF][START_REF] Sterman-Cohen | Rayleigh-Taylor instability in thin liquid films subjected to harmonic vibration[END_REF] . In the latter case, the amplitude of the vibration needs to be increased with the surface size. The maximum amplitude is set by the triggering of another instability called the Faraday instability which tends to destabilize fluid surfaces above a certain acceleration threshold [START_REF] Faraday | On a Peculiar Class of Acoustical Figures; and on Certain Forms Assumed by Groups of Particles upon Vibrating Elastic Surfaces[END_REF][START_REF] Douady | Experimental study of the Faraday instability[END_REF] . However, this threshold can be raised by increasing the fluid viscosity [START_REF] Kumar | Parametric instability of the interface between two fluids[END_REF] . Hence, the upside-down liquid volume can be large provided the viscosity is properly chosen.

The vertical vibration of a fluid also induces air bubbles to sink below a certain depth in the liquid defying the well-known Archimedes' principle [START_REF] Baird | Resonant bubbles in a vertically vibrating liquid column[END_REF][START_REF] Jameson | The motion of a bubble in a vertically oscillating viscous liquid[END_REF][START_REF] Sorokin | Motion of a gas bubble in fluid under vibration[END_REF][START_REF] Blekhman | Anomalous" phenomena in fluid under the action of vibration[END_REF][START_REF] Blekhman | Surface and volumetric effects in a fluid subjected to high-frequency vibration[END_REF][START_REF] Zen'kovskaja | Influence of vertical oscillations on bilaminar system with a non-rigid interface[END_REF][START_REF] Chelomei | Mechanical paradoxes caused by vibrations[END_REF] . This effect has been studied for industrial applications in gas holdup and mixing in bubble column reactors [START_REF] Elbing | Review of Bubble Column Reactors with Vibration[END_REF] .

Here, we investigate the effect of the vertical vibrations on the buoyancy of bodies immersed in levitating liquid layers and in particular at their lower interface. Our experimental setup consists of a Plexiglas container fixed on a shaker vibrated vertically at frequency ω/2π with an amplitude A (Fig. 1a). The container is filled with silicon oil or glycerol with high viscosity (typically ranging from 0.2 Pa.s to 1 Pa.s) to increase the Faraday instability threshold [START_REF] Wolf | Dynamic Stabilization of the Interchange Instability of a Liquid-Gas Interface[END_REF] . Though they have different physical properties, both liquids exhibit similar behavior provided that their viscosity is large enough. In particular, the wetting conditions appear to have limited influence thanks to dynamical effects on the contact line. Air bubbles are observed to sink when placed below a critical depth. This behavior which defies standard buoyancy can be explained by a simple model taking into account the kinetic force, also called the Bjerknes force [START_REF] Bjerknes | Fields of Force: Supplementary Lectures, Applications to Meteorology[END_REF] , exerted on the bubble in the oscillating bath [START_REF] Baird | Resonant bubbles in a vertically vibrating liquid column[END_REF][START_REF] Sorokin | Motion of a gas bubble in fluid under vibration[END_REF] (see Supplementary Materials and Movie 1 for details). By expanding an already sunken bubble, we create an air layer trapped below a levitated liquid layer (Fig. 1b and Supplementary Movie 2). Its lower interface is stabilized by the vertical shaking, preventing the release of the trapped air. This air layer is acting as a vertical spring loaded with the liquid mass placed upon it and driven by the shaker (Fig. 1c). It can be modeled by a driven damped harmonic oscillator 𝑧 +

2𝛤𝜔 &'( 𝑧 + 𝜔 &'(

) 𝑧 = 𝐴𝜔 ) cos(𝜔𝑡) with ω res the resonance frequency of the air layer and Γ the damping ratio due to the shearing induced by the relative motion between the levitating liquid layer and the bath walls (see Supplementary Information). In the laboratory frame, the normalized oscillation amplitude A l (ω)/A and its associated relative phase ϕ l (ω)-ϕ compared to the shaker clearly show the expected resonance behavior (Fig 1c and1d). The air layer thus enables the enhancement of the excitation amplitude of the shaker by more than one order of magnitude. Near resonance, the amplitude is high enough to excite the Faraday instability on both sides of the fluid layer (see inset Fig. 1d and Supplementary Movie 3). This resulting "rain" emitted from the lower interface induces a thinning of the fluid layer which can be avoided by simply reducing the excitation amplitude. Provided the spring-mass oscillation is properly tuned, there is no restriction in the number of sustained levitating layers which can be piled on top of one another (see Fig. 2e and Supplementary Movie 4).

As mentioned before, the vertical vibrations have a stabilizing effect on the lower fluid interface.

This can be interpreted as a Kapitza effect which consists of a dynamical stabilization of an inverted pendulum by vertical shaking [START_REF] Kapitza | Dynamic stability of a pendulum when its point of suspension vibrates[END_REF][START_REF] Krieger | Interfacial fluid instabilities and Kapitsa pendula[END_REF] . Solving the Bernoulli equation for the fluid shows that the interface height 𝜁 𝑘 at the spatial wavenumber 𝑘 behaves like an inverted pendulum. The spatial mode satisfies 𝜁 + 𝜔 0 (𝑘) 2 + 𝐴 𝑙 2 𝑘 2 2 𝜔 2 𝜁 = 0 with 𝜔 6 𝑘 ) = -𝑔𝑘 + (𝛾𝑘 ; )/𝜌 > the gravitocapillary dispersion relation with inverted gravity. Without vibrations, the oscillator is unstable for small enough k (ω 0 (k) 2 <0) leading to the Rayleigh-Taylor instability while large wave numbers are stabilized by capillarity. The last term in the equation arises from the modulation of the effective gravity. In gravitational regime, the stabilization is reached for wavenumbers satisfying 𝑘 > 2𝑔/𝐴 > ) 𝜔 ) (see Supplementary Information). The limited size L for the bath sets a maximum limit to the observed excitable wavenumber k>2π/L (only anti-symmetric modes satisfying volume conservation are considered). As a consequence, the stability of the interface is obtained for oscillating liquid velocities 𝑣 > = 𝐴 > 𝜔 above a critical velocity 𝑣 > * = 𝑔𝐿/𝜋 < 𝑣 E .

There seems to be no size limit for stabilization. The maximum levitated mass was 0.5 L in a 12x12 cm 2 container, and the max length achieved was 20 cm. The limitations in mass are only due to the shaker. In addition, no decay in time was observed and the layers remained stable for arbitrary long times. Figure 1f shows the oscillating velocity 𝑣 > * needed to stabilize baths with lengths L up to 18 cm (insets show views of levitating layers for L=2 cm and L=18 cm, see Supplementary Movie 5). Other effects might have an influence on the stability, such as friction or the flows formed at the boundary, but their influence is limited (see Supplementary Information).

We now focus on objects floating at the inverted interface of the levitating fluid layer. Archimedes' principle states that the upward buoyant force exerted on an immersed body, whether fully or partially submerged, is equal to the weight of the displaced fluid. Although this may seem counterintuitive, the transpose symmetric position at the lower interface (see Fig 2a) also exhibits an upward buoyant force equal to the weight of displaced liquid. Figure 2b shows the typical potential exerted on a floating body without taking into account the dynamic effects (see Supplementary information for details). The two equilibrium positions associated with each interface are clearly visible. However, while the upper position is stable, the lower is not: pushing the body down/up would make it fall/float to the upper interface. Taking into account the dynamical effect, i.e. the time averaged effect of the oscillations, provides an additional stabilizing dynamical potential around the two equilibrium positions (see inset Fig. 2b). Averaged small displacements Z b of the floater around the two equilibrium positions satisfies the same dynamical equation 𝑍 G + 𝜔 G ) 1 + 𝛼 𝑍 G = 0, ω b being the angular frequency associated with the buoyancy force and 𝛼 =

J K L ) M N J O )
the correction induced by the averaged dynamical effects (see Supplementary Information for details). While, the dynamical effects increase the stability of the equilibrium at the upper interface 𝜔 G ) > 0, 𝛼 > 0 , the unstable static equilibrium at the lower interface (𝜔 G ) < 0) is stabilized by the dynamical effects (α<0). It is interesting to note that similar dynamical stabilizations were observed with a washer mounted on a vibrated inverted pendulum [START_REF] Chelomei | Mechanical paradoxes caused by vibrations[END_REF] .

The floater stability is reached for liquid velocity 𝑣 > above a critical value given by 𝑣 G * = √2𝑔/|𝜔 G | < 𝑣 > . It is thus possible to have floating bodies with varying density above and below the levitating liquid layers (see Figure 2c). Hence the vibration not only gives stability of the lower horizontal interface of a liquid but also permits a vertical stabilization of the unstable equilibrium position that a floater would experience on such interface. This dynamical "anti-gravity" enables boats to float on both interfaces (Fig. 1e, see Supplementary Movie 6). 33 Note that the drag force induced by secondary flows due to recirculation in the liquid layer should not significantly change the vertical equilibrium position (see Supplementary Movie 2 and Supplementary Information).

The previous stability conditions suggest that the critical fluid velocity to stabilize a floater 𝑣 G * and the liquid layer 𝑣 > * are different and that sufficiently massive floaters should fall before the layer collapses (𝑣 G * > 𝑣 > * ). We performed the experiments in silicon oil with spherical floaters of increasing mass, the heaviest one being almost neutrally buoyant (Figure 3). Contrary to light floaters which fall with the liquid layer as the excitation amplitude is decreased, heavier floaters fall before the levitating layer gets destabilized (see inset Figure 3 and Supplementary Movie 7).

The theoretical critical velocity 𝑣 G can be exactly computed for spheres without any adjustable parameters (solid blue line). The expected range of masses for which 𝑣 G * > 𝑣 > * is consistent with the experimental findings (blue area in Fig. 3) and the values 𝑣 G * are in reasonable agreement.

Discrepancy occurs for almost neutrally buoyant floaters (as also for buoyant equilibrium position). In this limit, new phenomena are observed such as a small relative motion of floaters with respect to the surrounding fluid layer which seems to play a significant role in the floater equilibrium.

This counter-intuitive upside-down buoyancy phenomenon suggests that the stabilization of Rayleigh-Taylor instability through vibrations can be considered not only in itself but also as a playground for new experiments in unexplored conditions. Under this perspective, many phenomena that occur at the interface between air and liquids could be investigated and reformulated in this new exotic configuration such as transport and segregation.

Figure 1 Figure 2

 12 Figure 1 Levitating liquid layer stabilized by Kapitza effect. a, Experimental setup composed of a Plexiglas container of various sizes (up to 20 cm in width) attached on a vertically oscillating shaker with amplitude A and frequency ω/2π. The liquid is either glycerol or silicon

Figure 3

 3 Figure 3 Floater and liquid layer stability. Critical velocities for the stability of the liquid layer 𝑣 > * and for the floater 𝑣 G * as a function of the floater mass m b . The experiments are performed in a 4x5 cm 2 container with silicon oil and a spherical floater of 2.5 cm with various masses. The model (blue line) based on dynamic stabilization uses the liquid density measured experimentally ρ=0.92 kg/L and no adjustable parameter. Above a certain floater mass 𝑣 G * > 𝑣 > * (blue area) the floater can fall while the liquid layer remains stable. The error bars correspond to extremal values over 3 measurements. The dashed red line is the mean value of the 𝑣 > * measurements. Inset image sequences of the experiment showing the layer falling at 𝑣 > * for 𝑚 T> = 4.8 g and showing the floater falling at 𝑣 G * for 𝑚 G = 6.6 g (see Supplementary Movie 7).

Acknowledgments:

The authors are grateful to Suzie Protière, Arnaud Lazarus, Sander Wildeman and the colleagues and students of "Projets Scientifiques en Equipes" for insightful discussions.

The authors thank the support of AXA research fund and the French National Research Agency LABEX WIFI (ANR-10-LABX-24).

Competing interests: Authors declare no competing interests.

Additional information: Supplementary Information and Movies are available for this paper.

Correspondence and requests for materials should be addressed to E. F.

Data availability:

All the datasets generated during the current study are available in the Supplementary Information.


Alessandro Montaldo

Luca Fronda

Ihsen Hedhli

Gabriele Moser
email: gabriele.moser@unige.it.

Sebastiano B Serpico

Josiane Zerubia
email: josiane.zerubia@inria.fr

A CAUSAL HIERARCHICAL MARKOV FRAMEWORK FOR THE CLASSIFICATION OF MULTIRESOLUTION AND MULTISENSOR REMOTE SENSING IMAGES

Keywords: Multiresolution and multisensor fusion, causality, hierarchical Markov random models, Markov chain, semantic image segmentation

In this paper, a multiscale Markov framework is proposed in order to address the problem of the classification of multiresolution and multisensor remotely sensed data. The proposed framework makes use of a quadtree to model the interactions across different spatial resolutions and a Markov model with respect to a generic total order relation to deal with contextual information at each scale in order to favor applicability to very high resolution imagery. The methodological properties of the proposed hierarchical framework are investigated. Firstly, we prove the causality of the overall proposed model, a particularly advantageous property in terms of computational cost of the inference. Secondly, we prove the expression of the marginal posterior mode criterion for inference on the proposed framework. Within this framework, a specific algorithm is formulated by defining, within each layer of the quadtree, a Markov chain model with respect to a pixel scan that combines both a zig-zag trajectory and a Hilbert space-filling curve. Data collected by distinct sensors at the same spatial resolution are fused through gradient boosted regression trees. The developed algorithm was experimentally validated with two very high resolution datasets including multispectral, panchromatic and radar satellite images. The experimental results confirm the effectiveness of the proposed algorithm as compared to previous techniques based on alternate approaches to multiresolution fusion.

INTRODUCTION

Thanks to the potential offered by current and forthcoming space missions, very high spatial resolution (VHR) optical (e.g., Pléiades, WorldView-3, SPOT-6/7) and synthetic aperture radar (SAR; e.g., COSMO-SkyMed Seconda Generazione, TerraSAR-X, RADARSAT-2) instruments are now available to serve many applications. From the viewpoint of image analysis, the joint exploitation of the resulting data requires novel methods that operate with images collected by multiple sensors on the same area at multiple spatial resolutions. This is highly promising because it allows to benefit from data associated with different physical natures, frequencies, polarizations, etc., and from the tradeoff between a synoptic view at coarser resolutions and the spatial detail of finer resolutions. In this paper, the semantic segmentation [START_REF] Wang | Deep learning in object recognition, detection, and segmentation[END_REF][START_REF] Arnab | Conditional Random Fields Meet Deep Neural Networks for Semantic Segmentation: Combining Probabilistic Graphical Models with Deep Learning for Structured Prediction[END_REF] or dense supervised classification of images that are both multiresolution and multisensor is addressed. A major challenge in this joint multisensor-multiresolution scenario is the combination of heterogeneous statistics of the input images and of the need to characterize spatial information associated with different resolutions [START_REF] Hedhli | Classification of Multisensor and Multiresolution Remote Sensing Images Through Hierarchical Markov Random Fields[END_REF]. Trivial well-known solutions mostly use resampling procedures and do not attempt to capture the multiresolution structure of the data explicitly.

Here, an approach to address multiresolution fusion is proposed based on hierarchical latent Markov modeling [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF]. Within the family of structured output learning methods, Markov models postulated on planar or multilayer graphs are known as flexible and powerful stochastic models for spatial and possibly multimodal information [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF]. A common shortcoming of generic Markov random fields (MRFs) is their non-causality that generally leads to possibly time-consuming inference algorithms. The resulting non-convex problems are addressed with stochastic optimizers such as simulated annealing [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF] or graph-theoretic algorithms [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. In opposition to the case of classical MRFs, for causal probabilistic 1D models (e.g., Markov chains) efficient sequential techniques (e.g., the Baum-Welch algorithm [START_REF] Baggenstoss | A modified Baum-Welch algorithm for hidden Markov models with multiple observation spaces[END_REF]) are widely used for inference. Indeed, an effective approach to multiresolution and multisensor fusion addressed with Markov modeling is to postulate Markovianity on a hierarchical graph structure (e.g., systems of quadtrees [START_REF] Hedhli | Classification of Multisensor and Multiresolution Remote Sensing Images Through Hierarchical Markov Random Fields[END_REF]). However, despite hierarchical MRFs on quadtrees capture relations among sites located at different scales through the use of a Markov chain, a common limitation is that they do not explicitly characterize spatial dependences within the layer at each resolution [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF].

To join the benefits of both hierarchical and planar MRFs, the classification of multiresolution and multisensor remote sensing data is formulated in this paper through a causal hierarchical Markov framework. The contribution of the paper is twofold. First, we extend the approach developed in (Montaldo et al., 2019a, Montaldo et al., 2019b), in which hierarchical MRFs on quadtrees have been combined with planar Markov meshes and with various pixelwise probabilistic models, to a broader framework in which each layer of the quadtree is associated with a Markov model with respect to an arbitrary total order relation. The resulting framework is highly general and includes the Markov meshes in (Montaldo et al., 2019a, Montaldo et al., 2019b) as a special case. Secondly, within this generalized framework, we develop a novel algorithm that is based on a Markov chain model [START_REF] Fjortoft | Unsupervised classification of radar images using hidden Markov chains and hidden Markov random fields[END_REF] on each layer of the quadtree. In this model, Markov chains are formulated both across the scales of the tree and with respect to a 1D scan of the pixel lattice of each layer. This joint strategy benefits from the spatial information within each layer and inherently supports multiresolution fusion. A case-specific scan is defined on each layer of the quadtree by integrating a simple zig-zag trajectory and a Hilbert space-filling curve [START_REF] Abend | Classification of Binary Random Patterns[END_REF]. Methodologically, we prove that the proposed generalized framework ensures causality for the whole probabilistic graphical model, thus favoring efficient inference. Then, we derive the analytical formulation of the marginal posterior mode (MPM) inference algorithm for the proposed framework and we customize it to the special case of the aforementioned Markov chain. MPM is especially advantageous for classification and segmentation in a multiresolution scenario [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF] and its formulation for causal hierarchical Markov models explicitly relies on the causality of the models. In the proposed framework, the images coming from the various sensors are inserted in a quadtree based on their resolutions. On one hand, when distinct sensors provide data at different resolutions, multisensor fusion is accomplished naturally thanks to the adopted hierarchical topology. On the other hand, when data sources from distinct sensors are available at the same resolution, multisensor fusion is addressed by incorporating a nonparametric ensemble method within the proposed algorithm. Similar to (Montaldo et al., 2019a, Montaldo et al., 2019b), gradient boosting regression trees (GBRT) [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF] are employed in this role.

METHODOLOGY

Causal hierarchical Markov models

Among the main current trends of probabilistic graphical models in remote sensing, multiscale and multiresolution models are of primary importance. The rationale is that complementary information is appreciated at different spatial scales. While fine-scale observations provide lots of geometrical detail but are spatially heterogeneous and noise-sensitive, coarsescale observations capture only large image regions and covers but with strong robustness to noise and outliers. Many multiscale or multiresolution graphical models have been introduced on the basis of hierarchical graphs [START_REF] Willsky | Multiresolution Markov models for signal and image processing[END_REF], such as quadtrees [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. With this topology, hierarchical MRF models typically involve Markovianity along the scale axis. However, a limitation of a hierarchical MRF on a quadtree is that, while it captures relations among data at different scales through a Markov chain, it generally does not characterize spatial dependencies within the layer at each scale [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. Here, a causal hierarchical framework that postulates Markovianity not only among pixels belonging to pixel lattices with different resolutions but also among pixels in the same lattice is defined. This framework describes both cross-layer (multiscale) and intra-layer (spatial) dependence relations.

2.1.1 Planar causal Markov models Let R ⊂ Z 2 be a planar rectangular pixel lattice and let us define a total order relation on R. This means that this relation fulfills the following properties (r, s, t ∈ R): Antisymmetry: r s, s r =⇒ r = s Transitivity: r s, s t =⇒ r t Connexity:

at least one holds between r s and s r.

We shall also write r s to indicate that r s and r = s. Accordingly, the set {r ∈ R : r s} formalizes the idea Figure 1. Quadtree structure and notations on the tree of the "past" of site s ∈ R. Let us consider also a neighborhood relation to be defined in R consistently with the order relation , so that, intuitively, it makes sense to speak of the "past neighbors" or "recent past" of each site. We write r s to indicate that r is one of such past neighbors of s. Formally, {r ∈ R : r s} {r ∈ R : r s}, i.e., the past neighbors of s ∈ R are all included in the past of s but form a proper (strict) subset of its entire past. Normally, this subset is meant to be small to favor low computational burden in the inference.

Let a discrete random variable xs be attached to each site s ∈ R and let X = {xs}s∈R be the corresponding random field. A Markovianity property restricted to the past of each site holds for X if [START_REF] Abend | Classification of Binary Random Patterns[END_REF][START_REF] Devijver | Hidden Markov mesh random field models in image analysis[END_REF][START_REF] Willsky | Multiresolution Markov models for signal and image processing[END_REF]:

P (xs|xr, r s) = P (xs|xr, r s) ∀s ∈ S. (1)
This notion of causal Markovianity on a planar lattice extends the one discussed in (Montaldo et al., 2019a, Montaldo et al., 2019b), in which the order relation was referred to a 2D mesh. We also recall that it can be proved for several planar causal Markov models (including Markov meshes of second and third order and Markov chains) that the following factorization holds [START_REF] Abend | Classification of Binary Random Patterns[END_REF]:

P (X ) = s∈R P (xs|xr, r s), (2)
up to properly defining the behavior of X at the image borders.

2.1.2 Hierarchical causal Markov random fields Let now {S 0 , S 1 , . . . , S L } be a collection of pixel lattices organized as a quadtree, as shown in Fig. 1. This hierarchical structure is arranged so that the width and the height of S -1 are twice smaller than those of S ( = 1, 2, . . . , L). If ∈ {1, 2, . . . , L}, each site s ∈ S in layer S has a well-defined parent site s -∈ S -1 . If ∈ {0, 1, . . . , L -1}, each site s ∈ S has in S +1 four well-defined children sites, which we collect in the set s + ⊂ S +1 . These definitions imply a hierarchy on the whole tree S = S 0 ∪S 1 ∪. . .∪S L from the root S 0 to the leaves S L [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. If again a discrete random variable xs is associated with each s ∈ S, then X = {xs}s∈S is a hierachical MRF on the quadtree if the following property holds for all layers but the root ( = 1, 2, . . . , L) [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]:

P (X |X -1 ,X -2 , . . . , X 0 ) = = P (X |X -1 ) = s∈S P (xs|x s -), (3)
where X = {xs} s∈S ( = 0, 1, . . . , L). The first equality in (3) indicates that Markovianity is valid across the layers of the quadtree. The second equality expresses that the transition probability P (X |X -1 ) from the ( -1)th to the th layer factorizes on a pixelwise basis. Accordingly, a hierarchical MRF on a quadtree describes the dependence along the scale axis but not the spatial-contextual dependence associated with the random field X that is supported on each layer S .

2.2

The proposed framework 2.2.1 Model assumptions In this paper, we formulate a general framework for the joint fusion of multisensor and multiresolution images in a supervised classification scenario. The key-idea of the proposed framework is to generalize the approaches in (Montaldo et al., 2019a, Montaldo et al., 2019b) by introducing a hierarchical Markov model with respect to a quadtree topology and to an arbitrary total order relation on each layer of the tree. Considering a set of well-registered images acquired by distinct VHR sensors at different spatial resolutions on the same area, each image is included in a corresponding layer of the quadtree. Hence, the resolutions of the input images should be in a power-of-2 relation. This is a restriction, but given the resolutions of current VHR sensors, it can be easily met up to minor resampling [START_REF] Dikshit | An empirical investigation of image resampling effects upon the spectral and textural supervised classification of a high spatial resolution multispectral image[END_REF][START_REF] Inglada | Analysis of artifacts in subpixel remote sensing image registration[END_REF] and appropriate antialiasing filtering [START_REF] Alparone | Remote Sensing Image Fusion[END_REF][START_REF] Mallat | A wavelet tour of signal processing -The sparse way[END_REF]. The hierarchical model naturally fits the requirements of multiresolution fusion. Furthermore, the relations , , and are assumed to be defined within each planar layer of the quadtree. From a graphical viewpoint, this implies that every site s ∈ S of every layer S but the root ( = 1, 2, . . . , L) is linked to one parent s -∈ S -1 located in the upper layer and to the past neighbors {r ∈ S : r s} located in its own layer. A pixel in the root S 0 has no parent and only the causal neighbors, if any.

With the same notations as above, xs is now meant as the class label of site s ∈ S and a finite set Ω of classes is assumed to be defined by a training map (xs ∈ Ω, s ∈ S). For each s ∈ S, let Cs be the vector collecting the labels of all nodes linked to s (the "context" of s), i.e., Cs = {s -} ∪ {r ∈ S : r s} if s is not in the root (s ∈ S \ S 0 ), and Cs = {r ∈ S 0 : r s} if s ∈ S 0 is in the root. As the input multisensor and multiresolution images are framed in the quadtree, each s ∈ S is also associated with a feature vector ys. Thus, a random field Y = {ys}s∈S of observations is defined. If a site s is not in the leaf layer (s ∈ S \ S L ), then we shall denote as Ds the vector collecting the feature vectors of all sites that descend from s ∈ S in the quadtree. If s in the leaf layer (s ∈ S L ), we set Ds = ys. The proposed framework is defined by the following assumptions:

1. X is Markovian across the scales as in (3) and the following proportionality and factorization hold on all layers except the root ( = 1, 2, . . . , L):

P (X |X -1 , X -2 , . . . , X 0 ) = P (X |X -1 ) ∝ ∝ s∈S P (xs|xr, r s)P (xs|x s -). (4)
2. On the root, X 0 satisfies the factorization in (2):

P (X 0 ) = s∈S 0 P (xs|xr, r s).
(5)

3. Conditional independence holds for the observations given the class labels:

P (Y|X ) = s∈S P (ys|xs) = L =0 s∈S P (ys|xs). (6)
4. The label of each site, given the labels of the sites linked to it, only depends on the observations of its descendants and not on those of the other sites:

P (xs|Cs, Y) = P (xs|Cs, Ds) ∀s ∈ S. (7)
5. Given the observation field, the parent and neighboring labels of each site are conditionally independent:

P (Cs|Y) = P (x s -|Y) r s P (xr|Y) ∀s ∈ S \ S 0 r s P (xr|Y) ∀s ∈ S 0 (8)
6. The parent and neighboring labels of each site, when conditioned to the label of that site, are independent on the observations of its descendants and mutually independent:

P (Cs|xs, Ds) = P (Cs|xs) = (9) = P (x s -|xs) r s P (xr|xs) ∀s ∈ S \ S 0 r s P (xr|xs) ∀s ∈ S 0
These assumptions play different roles within the proposed framework. On one hand, methodologically, the framework is defined by Assumptions 1, 2, and 3. Assumption 1 expresses the key idea of the proposed combination of a hierarchical MRF on a quadtree and a causal planar Markov model with respect to a generic order relation . It implies that Markovianity holds across the scales and that the corresponding transition probabilities factorize so that the labels in each layer convey both contextual and multiscale dependencies. Accordingly, the model can characterize both multiresolution and spatial information. From this perspective, Assumption 2 is necessary for completeness, as it ensures that the analogous spatial factorization holds on the root as well. The conditional independence in Assumption 3 is a commonly accepted statement in the application of MRFs to image classification and is used in the proposed framework to favor tractability. On the other hand, Assumptions 4, 5, and 6 are not meant to formalize specifically desired behaviors of the proposed model. As we shall argue in more detail in the next section, they are basically technical conditions that favor the analytical tractability of the inference. They are similar in nature to conditional independence hypotheses that are often stated for mathematical convenience when operating with either hierarchical [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF] or planar MRFs [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF].

Regarding inference within the proposed framework, we use the MPM criterion, i.e., s ∈ S is assigned to the class ω ∈ Ω that maximizes P {xs = ω|Y} [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF]. MPM is especially appropriate for hierarchical MRFs because it penalizes errors according to the scale on which they are made, which is desirable to avoid error accumulation along the layers [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF]. We shall prove in the next section that the proposed framework is causal both spatially and across scales, which makes it possible to compute P (xs|Y) on each s ∈ S recursively and efficiently. As data associated with different resolutions are associated with the observations in distinct layers of the quadtree, this recursive inference inherently accomplishes their fusion, thus jointly classifying multisensor data supported at different resolutions. If two or more considered input sensors correspond to the same resolution, then their data contribute to the feature vectors in the same layer (stacked vector). In this case, multisensor fusion is obtained by integrating the GBRT ensemble of decision trees [START_REF] Friedman | Greedy function approximation: A gradient boosting machine[END_REF] in the MPM scheme (see Section 2.3). The flexibility of tree ensembles and their applicability to highly heterogeneous input feature vectors allows the proposed framework to address multisensor fusion in this case as well.

Methodological properties and MPM inference

Here, we prove the causality of the proposed hierarchical Markov framework and we derive a set of relations that allows computing the MPM functional P (xs|Y) of each site s ∈ S.

Operatively, the causality is expressed in terms of the factorization of the joint distribution of all observations and labels in terms of causal transition probabilities. Theorem 1. Under Assumptions 1-3, the joint distribution P (X , Y) of all class labels and feature vectors in the quad-tree is entirely defined by the parent-child transition probabilities {P (xs|x s -)} s∈S\S 0 , the past neighbor transition probabilities {P (xs|xr, r s)}s∈S, and the pixelwise data conditional likelihoods {P (ys|xs)}s∈S.

Proof. According to (3), Assumption 1 implies that:

P (X , Y) = P (Y|X )P (X ) = P (Y|X )P (X L , X L-1 , . . . , X 0 ) = P (Y|X )P (X L |X L-1 ) . . . P (X 1 |X 0 )P (X 0 ).

(10)

Plugging the factorizations (6) for P (Y|X ), (5) for P (X 0 ), and (4) for P (X |X -1 ) ( = 1, 2, . . . , L) into (10) yields:

P (X , Y) ∝ L =1 s∈S P (xs|xr, r s)P (xs|x s -)P (ys|xs)• • s∈S 0 P (xs|xr, r s)P (ys|xs). (11)
Therefore, P (X , Y) is entirely defined by the parent-child transition probabilities P (xs|x s -), the past neighbor transition probabilities P (xs|xr, r s) and the data conditional likelihoods P (ys|xs).

Theorem 2. Under Assumptions 1-6, for each site s ∈ S \ S 0 not in the root 1 :

P (xs) = x s - P (xs|x s -)P (x s -), (12)
Proof. First, ( 12) is a straightforward application of the total probability theorem, and the proof of ( 17) is the same as that reported in [START_REF] Laferté | Discrete Markov image modeling and inference on the quadtree[END_REF] for a hierarchical MRF. Then, the total probability theorem implies that:

P (xs|Y) = Cs P (xs|Cs, Y)P (Cs|Y), (18)
and ( 14) and ( 16) follow from plugging ( 7) and ( 8) into this equation. Moreover, thanks to Bayes theorem:

P (xs|Cs, Ds) ∝ P (Cs|xs, Ds)P (xs|Ds), (19)
where the proportionality constant does not depend on xs. If s ∈ S \ S 0 and if we plug (9) into this equality, we obtain:

P (xs|Cs, Ds) ∝ P (xs|Ds)P (x s -|xs) r s P (xr|xs), (20)
from which (13) follows due to Bayes theorem. The case of (15) is analogous.

Theorem 1 indicates that (X , Y) is a Markov chain with respect to the causal structure described above, which determines the causality of the proposed hierarchical Markov model. Hence, the proposed hierarchical Markov framework is causal both spatially and across scales, which allows an efficient recursive algorithm to be formulated for the MPM criterion. Based on Theorem 2, this criterion can be formulated in terms of three recursive steps. First, the prior probabilities P (xs) are initialized on the root and ( 12) is used to calculate them for each class on all sites of the quadtree through a top-down pass from the root to the leaves. Second, ( 17), ( 13), and (15) are used to compute P (xs|Cs, Ds) through a bottom-up pass from the leaves to the root and within the root. Third, ( 14) and ( 16) are used to derive P (xs|Y) through a second top-down pass within the root and from the root to the leaves.

Proposed Markov chain-based method

Different choices for the definition of the total order relation , meaning different choices of the set of pixels which represent the past of a site s ∈ S, lead to different algorithmic formulations. In (Montaldo et al., 2019a, Montaldo et al., 2019b), a Markov mesh model was considered, which, in the generalized formulation of the present paper, is equivalent to saying that r s if pixel r is a neighbor of s and is reached before s as the image is raster-scanned. Here, a Markov chain formulation is used, which implies that the past of a pixel is determined on a 1D basis according to a certain way of scanning all the pixels in a lattice. On each layer S of the quadtree ( = 0, 1, . . . , L), let us consider a scan trajectory, i.e., a sequence that visits every pixel once and moves from one pixel to one of its neighbors. Given a pixel s ∈ S , let us indicate as s ∈ S the previous pixel in the scan. We focus here on the special case in which r s if and only if r = s . Accordingly, ( 4) and ( 5) characterize the spatial behavior of X as a first-order Markov chain. In particular, ( 13) and ( 14) simplify as follows (s ∈ S \ S 0 ): 15) and ( 16) simplify analogously.

P (xs|x s -,
Two well-known ways to perform such a pixel scan include a zig-zag trajectory and the Hilbert space-filling curve. In the former case, a zig-zag scan over the lattice S is performed as shown in Fig. 2. With the exception of the pixels at the border of the lattice, s is the pixel adjacent to s and placed diagonally with respect to s in the direction of a zig-zag. Along the image borders, this trajectory aligns with the borders themselves. The latter scan option is based on the Hilbert space-filling curve over a power-two sized lattice S (see Fig. 3). The rationale for the introduction of this order relation is twofold. First, the Hilbert curve is useful because it gives a mapping between 1D and 2D spaces that preserves locality, a precious property for contextual classification. This means that two points which are close to each other in the 1D scan are also close to each other after folding. The converse may not always be true. Second, the Hilbert curve is particularly interesting in the case of a multiresolution data set as the recursive construction of the curve matches the power-of-2 relation of the layer sizes of the quadtree.

A typical drawback of causal Markov approaches is the possible presence of anisotropic corner-dependent artifacts. To prevent them, symmetrization procedures that ensure that, within the overall execution of the method, pixels are visited in a balanced and isotropic manner, are often necessary. In this paper, this drawback is addressed by defining the scan according to a symmetrized combination of the aforementioned zig-zag and Hilbert curves. The proposed scan is shown in Fig. 4 and results as a combination of four zig-zag scans and two Hilbert scans. It is accomplished by moving onto two different directions of the Hilbert curve and onto the zig-zag curves along the two diagonals, in both possible directions. Through this scan, each site is visited multiple times in a symmetric manner, which prevents the risk of geometrical artifacts that may occur while integrating contextual information in a 1D scan on a 2D pixel lattice. Indeed, if asymmetric scanning was used, directional artifacts [START_REF] Yousefi | Generating symmetric causal Markov random fields[END_REF].

The same approach used in (Montaldo et al., 2019a, Montaldo et al., 2019b) is used to define the transition probability P (xs|x s -) across consecutive scales. It consists of a parametric stationary model in which P {xs = ω|x s -= ω} = θ for all ω ∈ Ω, where θ is a hyperparameter of the method, and P {xs = ω|x s -= ω } is constant over all ω = ω (ω, ω ∈ Ω). Following the results obtained in [START_REF] Hedhli | A New Cascade Model for the Hierarchical Joint Classification of Multitemporal and Multiresolution Remote Sensing Data[END_REF] in the case of a multitemporal hierarchical MRF model, the impact of θ on the results of the method are expected to be minor. An analogous model is used for the transition probability P (xs|xs ) along the scan trajectory.

GBRT is used to estimate the pixelwise posteriors P (ys|xs) from the training samples of the classes. GBRT is an ensemble of decision trees that uses boosting [START_REF] Schapire | The strength of weak learnability[END_REF] to iteratively and adaptively manipulate the training set in order to generate multiple decision trees. A weight is assigned to each training sample, and at each iteration the weighted error on the training set is minimized by a decision tree. Then the weights are updated so that larger weight is given to wrongly classified samples and smaller weight is given to correctly classified samples. The final predictor is constructed by a weighted vote of the individual decision trees where the weights are given by the accuracy of the corresponding decision tree.

EXPERIMENTAL RESULTS

Data sets and experimental setup

The proposed method was experimentally validated with two VHR remote sensing datasets. The first one was collected over Port-au-Prince, Haiti, shortly after the 2010 earthquake. It consisted of an RGB GeoEye-1 image at 1.25 m resolution (1024 × 1024 pixels), an RGB-NIR (near infrared) Quick-Bird image at 2.5 m resolution, and a SAR COSMO-SkyMed stripmap image at 5 m resolution. For the COSMO-SkyMed stripmap modality, 5 m is the spatial resolution corresponding to 4-look data. The resolution of the multispectral channels of QuickBird is 2.4 m, so minor resampling was applied. In the case of GeoEye-1, the resolutions of the multispectral and panchromatic channels are 1.64 and 0.41 m, respectively. The 1.25 m image was obtained by resampling after pansharpening. Pansharpened data were used as this was the format in which GeoEye-1 imagery was provided by the Google Crisis Response initiative 2 . In all resampling processes, proper antialiasing filtering was applied [START_REF] Alparone | Remote Sensing Image Fusion[END_REF][START_REF] Mallat | A wavelet tour of signal processing -The sparse way[END_REF]. The thematic classes in the data set included "buildings," "asphalt," "containers," "vegetation," and "water." In the desired output map, these classes were meant to be discriminated at the finest observed resolution, i.e., 1.25 m. Yet, the proposed method also mapped them at the other observed scales, i.e., 2.5 and 5 m.

The second data set was acquired by the IKONOS optical multiresolution sensor in 2004 over the area of Alessandria, Italy. It consisted of a single-channel panchromatic image at 1 m resolution (1024×1024 pixels) and a 4-channel RGB-NIR multispectral image at 4 m resolution. The thematic classes in this data set were "urban," "agricultural," "rangeland," "forest," "bare soil," "wet soil," and "water." As in the case of the first data set, they were meant to be discriminated on the 1 m lattice, although the method also provided a classification map at 4 m as a by-2 http://www.google.com/relief/haitiearthquake/geoeye.html product. The urban area was meant to be discriminated as a whole, without separating possible subclasses associated with ground materials.

For both data sets, non-overlapping training and test sets were manually annotated in homogeneous areas by a specialist. No borders were taken into account in the training and test sets to avoid mixed pixels in the ground truth.

For comparison purposes, three previous approaches to multiresolution and multisensor fusion were considered. The first one consists in a planar MRF classifier in which unary potentials were based on the pixelwise predictions obtained by random forest [START_REF] Breiman | Random forests[END_REF] after resampling all the images to the finest resolution. The pairwise potential was the Potts model, and the energy minimization was through sequential tree-reweighted message passing (TRW-S) [START_REF] Kolmogorov | Convergent Tree-Reweighted Message Passing for Energy Minimization[END_REF]. The second one is the method in [START_REF] Moser | Multiresolution Supervised Classification of Panchromatic and Multispectral Images by Markov Random Fields and Graph Cuts[END_REF] for the classification of multiresolution optical images through MRF, graph cuts, and a linear mixture model. This method makes use of a Gaussian class-conditional model. In the application to SAR data, it was adapted by using a log-normal class-conditional model and by extending the technique accordingly. The third benchmark technique is the recently developed approach in (Montaldo et al., 2019a, Montaldo et al., 2019b) in which decision tree ensembles are combined with a causal hierarchical Markov model in which a Markov mesh is postulated at each resolution scale. The framework proposed in the present paper extends this previous technique. The algorithm developed here within the proposed framework differs from that in (Montaldo et al., 2019a, Montaldo et al., 2019b) especially in the modeling of spatial-contextual information within each layer of the quadtree (through Markov chains in the developed method and Markov meshes in (Montaldo et al., 2019a, Montaldo et al., 2019b)).

Experimental results and comparisons

The overall and class-wise accuracies obtained on the test sets are presented in Tables 1 and2. Spatial details of the classification maps are shown in Figs. 5 and6. The proposed technique obtained high accuracies on the test sets of both data sets, a result that suggests the effectiveness of the proposed causal hierarchical Markov framework and of the Markov chain-based algorithm developed within this framework. Remarkably, the accuracies obtained by the technique in (Montaldo et al., 2019a, Montaldo et al., 2019b) were very similar (identical in some cases) to those of the proposed technique. This similarity is consistent with the fact that both algorithms are special cases of the general framework formalized in the present paper and share the same multiresolution hierarchical structure based on a quadtree. On one hand, this similarity suggests that the spatial-contextual models based on Markov meshes and Markov chains, when integrated in the proposed framework and at least in the case of the considered data sets, allows for analogous capabilities in terms of class discrimination. From this perspective, large accuracy differences between the two approaches are generally not expected. On the other hand, this similarity in accuracy is contrasted by significantly different computational burdens. In the proposed chain-based algorithm, a pixel s is linked to the preceding pixel in the 1D scanning curve. Hence, ( 13) and ( 14), in their special case (21), involve looping over the three labels xs, x s -, and xs , which is computationally O(|Ω| 3 ). In the case of the previ-ous method in (Montaldo et al., 2019a, Montaldo et al., 2019b), the mesh formulation involves additional neighboring pixels, which leads, for instance, to O(|Ω| 4 ) or O(|Ω| 5 ) in the case of meshes with second-or third-order, respectively. Accordingly, the algorithm developed here exhibits significant computational advantages over the previous formulation in (Montaldo et al., 2019a, Montaldo et al., 2019b).

The other two previous approaches also obtained high classification performances, yet with lower overall and/or class-wise accuracies (e.g., for "forest" in the "Alessandria" data set). This further suggests the effectiveness of the proposed framework based on a hierarchical causal probabilistic graphical model as compared not only to the use of a planar MRF together with a traditional resampling process but also to a previous recent Markov model in which multiresolution fusion is formalized through linear mixtures.

A visual analysis of the classification maps confirm these comments (Figs. 5 and6). The map generated by the proposed algorithm was characterized by remarkable visual regularity without exhibiting oversmoothing issues. This is interpreted as a consequence of the contextual modeling obtained through Markovianity across both the scale and spatial domains.

CONCLUSIONS

A general hierarchical causal Markov framework has been proposed for the challenging problem of jointly classifying input image data that are both multiresolution and multisensor.

The proposed framework generalizes previous formulations that made use of Markov meshes by combining a multiresolution quadtree structure with a spatial Markovianity condition formulated with respect to a generic total order relation. The causality of the overall framework and the formulation of MPM have been analytically proven, thus generalizing the previous methodological results we proved in (Montaldo et al., 2019a). Within this framework, a novel classification algorithm based on a Markov chain model within each layer of the quadtree has been developed. Experimental results obtained with VHR satellite images collected by panchromatic, multispectral, and SAR sensors pointed out the capability of the proposed algorithm to achieve high class-wise and overall accuracy. This behavior is shared by the previous technique in (Montaldo et al., 2019a, Montaldo et al., 2019b), which can be seen as a special case of the general framework developed here. On one hand, these results confirm the capability of this framework to address multiresolution and multisensor fusion through a causal hierarchical Markov approach. A further desirable property of the proposed framework is that it allows classification maps to be obtained at all spatial resolutions in the input data set. On the other hand, the specific algorithm developed in the present paper exhibits significantly lower computational complexity than the previous one in (Montaldo et al., 2019a, Montaldo et al., 2019b). Indeed, the approach of associating each layer of the quadtree with a 1D Markov chain along Hilbert and zig-zag curves demonstrated advantageous from a computational perspective thanks to its lower complexity and correspondingly shorter execution times. This reduction in the order of complexity allows to address more efficiently larger images and allows especially to consider a larger number of classes. Higher accuracies were also obtained by the proposed algorithm as compared both to a traditional resampling-based approach and to an advanced multiresolution model based on linear mixture [START_REF] Moser | Multiresolution Supervised Classification of Panchromatic and Multispectral Images by Markov Random Fields and Graph Cuts[END_REF]. Furthermore, these previous approaches allow a classification map to be generated at the finest spatial resolution in the input data set but not at the other resolutions involved (unless case-specific resampling procedures are applied).

Future extensions of this work may involve integrating its multiresolution structure with convolutional neural networks, whose convolution and pooling layers intrinsically extract multiscale data representations [START_REF] Goodfellow | Deep learning[END_REF], and with adaptive multiresolution topologies based for instance on superpixel segmentations at different scales [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF].

Figure 2 .

 2 Figure 2. Zig-zag scan

  xs , Ds) ∝ P (xs|Ds)P (xs) -2 • (21) • P (xs|x s -)P (x s -)P (xs|xs )P (xs ), P (xs|Y) = x s -,x s P (xs|x s -, xs , Ds)P (x s -|Y)P (xs |Y), and (

Figure 3 .

 3 Figure 3. Hilbert space-filling curve

  Steps 1, 2, and 3 of the scan Steps 4, 5, and 6 of the scan

Figure 4 .

 4 Figure 4. Symmetrized scan introduced. It consists of a combination of four zig-zag scans and two Hilbert scans (6 steps illustrated in the figure).

  Figure 5. "Haiti" data set. Details of (a) the optical image at 1.25-m resolution and of the classification maps obtained using (b) a single-resolution MRF after resampling, (c) an adaptation of the method in (Moser et al., 2016) to the case of input multisensor optical-SAR imagery, and (d) the proposed algorithm. Class legend: containers, vegetation, asphalt, buildings, water.

  Figure 6. "Alessandria" data set. Details of (a) the IKONOS panchromatic image at 1-m resolution and of the classification maps obtained using (b) a single-resolution MRF after resampling, (c) the method in (Moser et al., 2016), and (d) the proposed algorithm. Class legend: urban, agricultural, rangeland, bare soil, forest, water, wet soil.

  For each root site s ∈ S 0 : For each site s ∈ S \ S L not in the leaf layer:

	P (xs|Cs, Ds) ∝	P (xs|Ds) P (xs) |Cs|	r s	P (xs|xr)P (xr),	(15)
	P (xs|Y) =	P (xs|Cs, Ds)		P (xr|Y).	(16)
	Cs				r s
	P (xs|Ds) ∝ P (xs|ys)	t∈s + x t	P (xt|Ds)P (xt|xs) P (xt)	.

P (xs|Cs, Ds) ∝ P (xs|Ds)P (xs|x s -)P (x s -)P (xs) -|Cs| • • r s P (xs|xr)P (xr), (13) P (xs|Y) = Cs P (xs|Cs, Ds)P (x s -|Y) r s P (xr|Y). (14)

1 We denote as |A| the cardinality (i.e., the number of elements) of a finite set A.

Table 1 .

 1 "Haiti" data set. Overall accuracy (OA) [%], Cohen's κ, and class-by-class accuracies [%] on the test set for the proposed algorithm and the compared benchmark methods.

	resolution 1.25 m	containers vegetation asphalt buildings	sea	OA	κ
	Single-res. MRF after resampling	63.97	74.81	98.64	99.30	76.24 94.58 0.8880
	Adaptation of (Moser et al., 2016)	75.06	85.32	95.89	97.75	99.36 95.11 0.9010
	Method in (Montaldo et al., 2019b)	87.08	33.27	95.17	99.04	97.18 96.36 0.9124
	Proposed method	86.97	34.57	97.82	99.36	100	96.90 0.9057
	resolution 2.5 m	containers vegetation asphalt buildings	sea	OA	κ
	Method in (Montaldo et al., 2019b)	86.31	32.71	94.87	100.00	96.04 96.63 0.9276
	Proposed method	86.92	33.08	97.29	100.00	100	97.00 0.9218
	resolution 5 m	containers vegetation asphalt buildings	sea	OA	κ
	Method in (Montaldo et al., 2019b)	87.98	25.39	96.55	100.00	88.88 96.01 0.9288
	Proposed method	87.98	25.39	96.55	100.00	88.88 96.01 0.9288

Table 2 .

 2 "Alessandria" data set. Overall accuracy (OA) [%], Cohen's κ, and class-by-class accuracies [%] on the test set for the proposed algorithm and the compared benchmark methods.

	resolution 1 m	urban	agri	rangeland forest water wet soil bare soil	OA	κ
	Single-res. MRF after resampling	98.58 99.12	92.23	36.98	100	98.30	96.82	96.03 0.9453
	Method in (Moser et al., 2016)	99.70 99.21	97.34	64.92	100	100	99.66	98.60 0.9804
	Method in (Montaldo et al., 2019b)	100	99.07	99.65	100	99.65	100	100	99.12 0.9976
	Proposed method	99.54 98.26	99.91	86.40 99.80	100	100	98.58 0.9853
	resolution 2 m	urban	agri	rangeland forest water wet soil bare soil	OA	κ
	Method in (Montaldo et al., 2019b) 99.54 98.06	99.91	78.74 99.53	100	100	98.01 0.9825
	Proposed method	99.49 98.08	99.91	79.52 99.84	100	100	98.03 0.9829
	resolution 4 m	urban	agri	rangeland forest water wet soil bare soil	OA	κ
	Method in (Montaldo et al., 2019b) 99.09 96.81	99.82	64.86 98.72	100	100	96.36 0.9697
	Proposed method	99.09 96.81	99.82	64.86 98.72	100	100	96.36 0.9697
	concepts

ACKNOWLEDGMENTS

The authors would like to thank the Italian Space Agency (ASI), DigitalGlobe Inc., and Google Crisis Response for providing the COSMO-SkyMed, QuickBird, and GeoEye-1 images used for experiments, respectively.


Bahram Jafrasteh
email: bahram.jafrasteh@inria.fr

Isabelle Manighetti

Josiane Zerubia
email: josiane.zerubia@inria.fr

GENERATIVE ADVERSARIAL NETWORKS AS A NOVEL APPROACH FOR TECTONIC FAULT AND FRACTURE EXTRACTION IN HIGH-RESOLUTION SATELLITE AND AIRBORNE OPTICAL IMAGES

Keywords: Remote sensing, Deep learning, Curvilinear feature extraction, Image processing, Generative adversarial networks, High resolution, Tectonic fault and fractures, Fault mapping

We develop a novel method based on Deep Convolutional Networks (DCN) to automate the identification and mapping of fracture and fault traces in optical images. The method employs two DCNs in a two players game: a first network, called Generator, learns to segment images to make them resembling the ground truth; a second network, called Discriminator, measures the differences between the ground truth image and each segmented image and sends its score feedback to the Generator; based on these scores, the Generator improves its segmentation progressively. As we condition both networks to the ground truth images, the method is called Conditional Generative Adversarial Network (CGAN). We propose a new loss function for both the Generator and the Discriminator networks, to improve their accuracy. Using two criteria and a manually annotated optical image, we compare the generalization performance of the proposed method to that of a classical DCN architecture, U-net. The comparison demonstrates the suitability of the proposed CGAN architecture. Further work is however needed to improve its efficiency.

INTRODUCTION

Fractures and faults are ubiquitous on Earth. Fractures are generally small-scale (lengths < 10 -3 -10 3 m), shallow (a few km depth) cracks in the rocks, whereas faults are larger features (lengths 10 -3 -10 3 km, widths 1 -10 2 km from ground surface to depth) accommodating significant cumulative displacements (10 -3 -10 3 km) over millions of years (e.g., Manighetti et al., 2001a). Faults have a 3D architecture, including fault elements of different lengths and a myriad of fractures off the principal fault plane (e.g., [START_REF] Perrin | Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth[END_REF]. When submitted to large stresses, as those imposed by Plate Tectonics at the global scale, faults may produce earthquakes. Because these earthquakes can be devastating, decades of scientific works have been done to understand them and provide clues to anticipate them. Among key information, the way fractures and faults are distributed in a region partly controls the location, distribution and magnitudes of the earthquakes in that region. Locating faults and fractures accurately is thus of upmost importance to assess seismic hazard (e.g., [START_REF] Benedetti | Earthquake synchrony and clustering on Fucino faults (Central Italy) as revealed from in situ 36Cl exposure dating[END_REF][START_REF] Walpersdorf | Present-day kinematics and fault slip rates in eastern Iran, derived from 11 years of GPS data[END_REF]. This can be done by identifying the traces that faults form as they intersect the ground surface. These surface traces have specific properties that can be recognized visually. Therefore, many studies have been conducted to identify fault and fracture surface traces in many regions worldwide. Trace identification was done visually, from direct observation in the field or from analysis of remote sensing image data, while mapping of the fault and fracture traces, commonly called "annotation", was done manually (e.g., [START_REF] Manighetti | Propagation of rifting along the Arabia-Somalia plate boundary: Into Afar[END_REF]Manighetti et al., , 2001b[START_REF] Manighetti | Recovering paleoearthquake slip record in a highly dynamic alluvial and tectonic region (Hope Fault, New Zealand) from airborne lidar[END_REF][START_REF] Manighetti | Repeated giant earthquakes on the Wairarapa fault, new Zealand, revealed by Lidar-based paleoseismology[END_REF].

However, manual annotation on image data requests a tremendous time, and relies on expert knowledge, which might not be always available. This calls for the need to develop a fast, reliable and accurate automatic fault identification and mapping * Corresponding author method in image data, that is derived from expert knowledge yet can be operated without further inputs of the experts. Such methods have been developed in earlier works [START_REF] Cartabia | Lineaments recognition for potential fields images using a learning algorithm for Boltzmann machines[END_REF][START_REF] Mavrantza | Implementation and evaluation of spatial filtering and edge detection techniques for lineament mapping: case study -Alevrada[END_REF][START_REF] Aghaee Rad | Machine learning of lineaments from magnetic, gravity and elevation maps[END_REF][START_REF] Farahbakhsh | Computer vision-based framework for extracting tectonic lineaments from optical remote sensing data[END_REF], and are generally divided in two categories: unsupervised and supervised. The latter takes advantage of the expert knowledge and is thus expected to generate reliable results; however, it requests annotations produced by the experts, what may be a long and difficult task. On the contrary, the unsupervised methods do not involve any annotation inputs from the experts, and are thus of more simple use; however they are more prone to mistakes such as detecting irrelevant features. In the domain of fault detection, the unsupervised methods might for instance mistakenly consider any linear feature, such as image edges or roads, as fault or fracture traces. Correcting these errors would then request significant post-processing by the experts.

Tectonic faults have curvilinear traces at the ground surface. We claim here that these curvilinear traces can be extracted from image data using a combination of computer vision and machine learning methods, provided that we have access to expert annotations. We are thus conducting supervised learning. Supervised learning is a machine learning task that learns how to model a series of input variables to one or more output variables. One of the most common non-linear methods of supervised learning is called Artificial Neural networks (ANN). ANN, inspired by biological neurons, is an information processing unit having artificial neurons and "hidden layers". Each hidden layer can contain one or more neurons. The neurons are dedicated to collapse and map the received input variables into one single variable. For that, they are conducting different operations to the input data such as attributing them some weights, modifying them through bias, or enhancing them through linear or non-linear functions called activation functions. The neuron output is then sent to one or more neurons of the next hidden layer. These operations are repeated until the total information reaches the last output layer. Weights and biases are the main variables that need to be adjusted during the calculations, and these adjustments are called the ANN training. This training is usually done with a back-propagation algorithm (details can be found in [START_REF] Leung | The complex backpropagation algorithm[END_REF]). An ANN with more than one hidden layer is known as a Deep Neural Network (DNN). Image processing techniques can be combined with DNN to improve the performance of both methods. One of the most commonly used image processing techniques is image convolution. Image convolution is a process through which one pixel is added to its local neighboring pixels, using a weighted matrix known as kernel or filter. If we replace the DNN neurons with an image convolution operation, then the weights of the network are replaced by convolution kernels. This architecture, which benefits from the high computational capacity of neural networks and the local feature extraction capability of image convolution, is called Deep Convolutional Network or DCN. Here we propose to use a particular architecture of DCN: a Generative Adversarial Network (GAN) [START_REF] Goodfellow | Generative adversarial nets[END_REF]. GAN, inspired from game theory, combines two DCNs in a min-max alternating fashion: A first network, called Generator (G), is trained to generate synthetic images, also called fake images in GAN, that resemble the annotated ones; then a second network, called Discriminator (D), is trained to discriminate, through comparisons between synthetic and annotated images, which synthetic images produce the most reliable results. Because we need to present an annotated ground truth image to the GAN system, the architecture is called a "Conditional generative adversarial network" or CGAN. The CGAN approach has been used in earlier works to successfully extract roads and buildings from aerial and spatial images [START_REF] Pan | Building extraction from high-resolution aerial imagery using a generative adversarial network with spatial and channel attention mechanisms[END_REF][START_REF] Zhang | Aerial image road extraction based on an improved generative adversarial network[END_REF], but, to our knowledge, it has never been applied to fault and fracture detection.

Therefore, in this study, we develop a procedure to adapt CGAN for fault and fracture extraction in optical image data. We also investigate the relative efficiency of CGAN and DCN for fault and fracture extraction.

AVAILABLE APPROACHES

As said earlier, the available approaches are classified in two main categories: unsupervised [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF][START_REF] Jeong | Stochastic model for curvilinear structure reconstruction using morphological profiles[END_REF][START_REF] Jeong | Progressive Tree-like Curvilinear Structure Reconstruction with Structured Ranking Learning and Graph Algorithm[END_REF] and supervised methods [START_REF] Parikh | Application of neural networks to pattern recognition problems in remote sensing and medical imagery[END_REF][START_REF] Cartabia | Lineaments recognition for potential fields images using a learning algorithm for Boltzmann machines[END_REF][START_REF] Aghaee Rad | Machine learning of lineaments from magnetic, gravity and elevation maps[END_REF].

Most unsupervised methods rely on edge detection. Edge detection aims to identify pixels in an image where intensity and brightness change sharply. Algorithms based on edge detection have been widely used to extract lineaments, including faults, from satellite images [START_REF] Parikh | Comparison of backpropagation neural networks and statistical techniques for analysis of geological features in landsat imagery[END_REF][START_REF] Mavrantza | Implementation and evaluation of spatial filtering and edge detection techniques for lineament mapping: case study -Alevrada[END_REF][START_REF] Marghany | Lineament mapping using multispectral remote sensing satellite data[END_REF][START_REF] Corgne | An integrated approach to hydro-geological lineament mapping of a semi-arid region of West Africa using Radarsat-1 and GIS[END_REF][START_REF] Rahnama | TecLines: A MATLABbased toolbox for tectonic lineament analysis from satellite images and DEMs, part 1: Line segment detection and extraction[END_REF]. Introducing pre-processing steps such as histogram equalization [START_REF] Marghany | Lineament mapping using multispectral remote sensing satellite data[END_REF] or directional filters [START_REF] Mavrantza | Implementation and evaluation of spatial filtering and edge detection techniques for lineament mapping: case study -Alevrada[END_REF][START_REF] Marghany | Lineament mapping using multispectral remote sensing satellite data[END_REF][START_REF] Soto-Pinto | A new code for automatic detection and analysis of the lineament patterns for geophysical and geological purposes (ADAL-GEO)[END_REF][START_REF] Borisova | Algorithms for lineaments detection in processing of multispectral images[END_REF][START_REF] Farahbakhsh | Computer vision-based framework for extracting tectonic lineaments from optical remote sensing data[END_REF] has helped to improve the lineament extraction. However, edge detection generally fails detecting faults accurately because it identifies as faults other irrelevant linear features. ANN has also been applied to lineament detection in optical images [START_REF] Parikh | Application of neural networks to pattern recognition problems in remote sensing and medical imagery[END_REF][START_REF] Cartabia | Lineaments recognition for potential fields images using a learning algorithm for Boltzmann machines[END_REF][START_REF] Aghaee Rad | Machine learning of lineaments from magnetic, gravity and elevation maps[END_REF]. The approach can handle multiple inputs and outputs and clearly provides more satisfying results. Among these studies, Aghaee Rad, 2019 used a DCN to compute the pixel probability in an image to be a fault. However, they did not use optical images.

In this study, we go beyond these prior works and propose a novel approach to accurately extract lineaments, here faults and fractures, from optical images, based on a generative adversarial neural network.

NEW METHODOLOGY

We have developed a conditional GAN structure (Figure 1). The network consists of two DCN sub-networks: a Generator network G and a Discriminator network D. The Generator, also called Segmentor, annotates each pixel of the input image with a generated label probability. The Discriminator, also called Classifier, compares the generated labels for each pixel in the image with the ground truth annotated image and sends its feedback to the Generator.

Deep Convolutional Neural networks

Generally, three datasets are introduced in a neural network architecture: a training dataset, made to learn the DCN with expert annotations of the elements of concern in a part of the input image; a test dataset made to quantify the ability of the trained DCN (generalization ability) to predict the elements of concern in another unseen part of the image; and a previously unseen validation dataset made to both check the generalization performance as training is progressing, and to stop the training before we over-train the DCN ("overfitting"); overfitting indeed happens when the network is trained so much on a set of images that it can only reproduces that images, with high accuracy. As the computational cost to feed a large size image is very high, the training is done on smaller sub-images extracted from the source image. Each sub-image is a matrix of pixel values that form the network input. We use here a DCN, called U-net (upper right part of Figure 1) [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. U-net has two major parts, a first one which is a contraction path called Encoder, dedicated to identify the elements of concern, and a second expansion part, called Decoder, dedicated to locate accurately these identified elements within the image. The first layer of the U-net includes the input image as a multitude of matrices. The next layer, which is the first hidden layer, convolves (i.e., applies a filter) each input matrix to highlight some basic features in the image (such as edges, sharp brightness changes, etc.). The convolution is done using a sliding window (called stride) of 1 pixel that rolls over the sub-image, and a filter weight, as shown in Figure (2). This process reduces the final size of the image (downsampling). In the next step, the convolved image (also called "feature") is presented to the next layer, and the same operation is performed with more complex filters. The process is repeated until the multi-convolved and downsampled features reach the end of the contraction path. From this point forward, the image sustains a transposed convolution, that is, similar to the convolution but in a backward direction (Figure ( 2)). The objective here is to find the accurate position of the identified elements in the image. It results that the image size is increased across the successive layers (upsampling), until it is back to its original size. The main objective of the training phase is to optimize the weights of the convolutions (Figure ( 2)). As these weights are optimized, the network generates images that resemble more closely to the annotated ground truth. This optimization is done to minimize, using a loss function, the difference between the ground truth and the network output labels. The most common optimization is done by back-propagation using a Stochastic Gradient Descend algorithm [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF]. Over the entire calculation, different procedures are applied to accelerate the training (such as mini-batch and batch normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF] ). Skip connections are used to transfer the information between the encoder and the decoder parts. They thus connect the spatial feature maps in the Encoder with the reconstructed feature maps in the Decoder, and improve the generalization ability of the U-net [START_REF] Ye | Understanding geometry of encoder-decoder CNNs[END_REF].

One of the commonly used activation functions in the U-net architecture is Leaky ReLu [START_REF] Xu | Empirical evaluation of rectified activations in convolutional network[END_REF]. Leaky ReLU maps negative value of x to a linear function with slope λ and positive values of x to f (x) = x. It is defined as :

f (x) = λx x < 0 x x ≥ 0 (1)
Here, instead of using one DCN architecture such as U-net, we use two DCNs organized as a generative adversarial network, similarly to a two players game. In this architecture, one of the networks, called Generator, learns to produce images as close as possible to the ground truth ones. The other network, called Discriminator, has access to both the ground truth images and the images produced by the Generator, and learns to discriminate them. The Generator changes its parameters according to the feedbacks it receives from the Discriminator, making it to progressively learn to produce images that resemble more and more to the ground truth ones.

GAN versus conditional GAN

Let's assume we have a dataset with N training images xn and the corresponding ground truth label yn, with probability distributions P datax and P datay , respectively. Then, the conventional loss function [START_REF] Goodfellow | Generative adversarial nets[END_REF] of the GAN network is defined as follows:

min D max G LGAN (D, G) = Ex∼P datax [logD(x)] + Ez∼P z [log(1 -D G(z) ] (2)
where z is an input noise variable drawn from a probability distribution Pz. Conditional version of GAN can be constructed by feeding extra information yn to both the Generator and Discriminator networks, instead of noise variable. One can define the following loss function for the Conditional GAN (CGAN):

min D max G LcGAN (D, G) = Ex∼P datax ,P datay [logD (x, y) ] + Ex∼P datax [log(1 -D (x, G(x)) ] (3)
The first term in the above equation maximizes the probability to correctly predict the output variable y. The second term is used to train the Generator network G to accurately segment the image and fool the Discriminator network D.

The two networks G and D are adversarially trained, that is, the Generator network minimizes the difference between the network-annotated images and the ground truth, while the Discriminator network separates the network-annotated images from the ground truth ones.

Conditional GAN for optical image segmentation

Let's H and W be the height and the width of the n th image xn, respectively. The dimension of the input image is H × W × 3, as we use three bands of the image, and the dimension of the output segmented image is

H × W × 1.
The proposed loss function with N training samples is defined as:

min D max G LcGAN (D, G) = 1 N N n=1 L l=1 ||f l D (xn, G(xn)) -f l D (xn, yn)|| + 1 N N n=1 ||G(xn) -yn|| 2 (4)
Where f l D represents the features extracted from the l th layer of the Discriminator network. The first part of the equation refers to the mean absolute error (MAE), called the L1 loss function and the second part represents the mean squared error (MSE), called the L2 loss function. [START_REF] Xue | Segan: Adversarial network with multi-scale l 1 loss for medical image segmentation[END_REF] mentioned that the first term of the above loss equation captures image features at different scales, thus hierarchical features. We defined the last term to empower the Generator network in extracting faults and fractures.

Network architecture and training specifications

As said earlier, the Generator network is a fully connected convolutional neural network with an encoder followed by a decoder. Transposed convolution layers are selected to do upsampling. The parameters of convolutional and transposed convolutional layers are as follows: kernel size: 4 × 4, stride: 2, activation function: leaky ReLU. The network takes advantage of the skip connection layers from U-net to recover the spatial information that can be lost otherwise from down-sampling at various resolutions [START_REF] Drozdzal | The importance of skip connections in biomedical image segmentation[END_REF]. The Discriminator network uses the encoder part of the Generator network and merges hierarchical features of the image to capture both long-range and short-range spatial relations among pixels (Figure 1). The parameters of this network are as follows: kernel size: 4 × 4, stride: 2, activation function: leaky ReLU. The back-propagation method is used to optimize the network parameters. There are two steps for the training: in the first step, we fix the Discriminator parameters and train the Generator. In the next step, we fix the parameters of the Generator and train the Discriminator using the same gradient as that obtained in the previous step. The two networks are competing to alternatively maximize and minimize the same loss function. As the epochs (i.e., number of iterations) of training increase, the two networks learn to better predict the output. In the final steps, the Generator produces probability maps that resemble the ground truth images manually annotated by the experts. between the ground truth image and the images produced by the Generator. Training specifications are as follows: batch size number: 64; λ = 0.2; optimization algorithm: Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF] with initial learning rate: 0.002, β1 = 0.5 and β2 = 0.99; number of epochs: 200. However, we allow the network to continue the training until the performance stops increasing over twenty consequent epochs.

1 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 1 0 1 0 1 1 0 0 0 1 0 1 1 0 0 0 1 0 0 1 1 0 0 0 1 1 1 0 0 0 0 1 1 0 1 0 1 0 1 1 0 0.

EXPERIMENTAL STUDIES

Dataset: We use a georeferenced optical four band image with spatial resolution 5 × 10 -4 m and size 20164 × 6807, constructed with photogrammetry from hundreds of hand-held camera pictures taken on the field (Figure 3(a)). Here, we use only the red (R), green (G) and blue (B) bands of the image. The image encompasses faults and fractures cutting a granite rock. The distorted or incomplete parts of the image at the site corners were not used. The entire image is divided in three regions, based on the geographical location, for training, validation and test, including 384, 94, 222 sub-images, respectively. This type of division which generates test sub-images that are far from the training ones is also recommended by [START_REF] Jafrasteh | Comparison of machine learning methods for copper ore grade estimation[END_REF]Bastani et al., 2018. Figure 3(c)-3(d) shows one of the patches in the image. Figure 3(b) shows the available expert annotations. As these annotations are not homogeneous (parts of the image best annotated than others), we imposed a threshold to the input images: an image is presented to the network only if the number of its annotated pixels is higher than one hundred. Furthermore, there are different uncertainties in manual fault mapping: one is that actual fault traces cannot be located with a greater precision than a few pixels, even by the best experts; a second one is that some fault traces are uncertain, even for the best experts (see different colors in Figure 3(d), the red lines represent major certain faults, yellow lines are for minor certain faults, and green color with dashed lines represent uncertain faults). To address the former issue, we applied a Gaussian filter with σ = 1.0 to the annotated fault traces. The value of σ was selected based on trial and error in range 0.5-2.5, to find the best compromise between actual line thickness and efficient training of the model. To address the second issue, we basically considered two pixel classes: certain fault, and not-a-fault Data Augmentation: Since the training samples are few, 384, we used different augmentation techniques to artificially increase the number of training samples, while making the algorithm invariant to changes in image brightness, contrast, fault orientations and positions, etc. [START_REF] Perez | The effectiveness of data augmentation in image classification using deep learning[END_REF]. More specifically, we applied to the input images some horizontal and vertical reversals, crops and rotations, along with brightness, contrast, hue, and saturation changes [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Taylor | Improving deep learning using generic data augmentation[END_REF]). An input image sustains each of these changes with a probability of 0.5.

Baselines:

We compare the performance of our proposed method with that obtained with classical U-net. The architec- ture of U-net and the network parameters are the same as in our Generator network (Figure 1).

Metrics:

We use a confusion matrix (Table 1) to classify the pixels as True Positive (TP), True Negative (TN), False Negative (FN) and False Positive (FP). TP and TN are predicted pixels in agreement with the same manual annotations from an expert (being a fault and not being a fault, respectively). FP and FN are predicted pixels in disagreement with the manual annotations (prediction of a fault and of absence of a fault where the manual annotations say the opposite). We investigate the generalization performance of our method using Dice similarity coefficient and Intersection Over Union (IOU). The higher the values of these two criteria the better is the accuracy of the proposed method.

Dice similarity coefficient [START_REF] Sudre | Generalised dice overlap as a deep learning loss function for highly unbalanced segmentations[END_REF] is:

Dice = 2T P 2T P + F N + F P . (5)
IOU [START_REF] Jaccard | The distribution of the flora in the alpine zone. 1[END_REF] is :

IOU = T P T P + F N + F P . (6)
A perfect match between ground truth and predicted maps would make these two criteria equal to one.

Ground truth

Predicted pixel value Figure 4 shows a patch of the optical image in the test region, and compares the expert annotations with our predictions. The latter are presented as probabilities. We see that our predictions well produce the major faults. They also detect many of the more minor faults, even though their continuity is not totally well predicted. Many short fractures are predicted at a scale that goes beyond the annotated map. They might identify real, unmapped features, but this needs to be examined in subsequent work by the experts. We also note that the network mistakenly identified the edges of a rock inclusion as a fault set. This problem arises from the training set not including such inclusion. Therefore this problem will be solved easily in subsequent work by introducing images containing rock inclusion into the training set. Altogether our results demonstrate that Conditional GAN architecture with careful specifications can be employed to analyze optical images and extract curvilinear fault traces accurately. To the best of our knowledge, this is the first time that GAN architecture is adapted to extract geological faults and fractures from optical images. The main advantage of this approach compared to unsupervised methods is its ability to capture both linear and non-linear structures in a meaningful manner.

RESULTS AND DISCUSSION

CONCLUSIONS

In this study, we proposed a novel methodology to extract fault and fracture traces in optical images. The methodology is based on conditional GAN architecture. We conditioned the network from a manually annotated image. As the actual fault traces have a certain "thickness", we represented them with a Gaussian filter. Then, we introduced these classified annotated images to CGAN. The CGAN takes advantage of a novel loss function to improve the performance and extract hierarchical features from the input image. Using the Dice and IOU criteria, we compared the performance of our CGAN approach with that of the most classical U-net architecture, both applied to a large size and high resolution optical image. The higher values of these criteria for the CGAN approach demonstrate its higher performance. More work still needs to be done however, especially to include more fault classes reproducing their natural hierarchy (major faults, more minor faults, faults uncertain at different degrees, etc.), and to take into account that manual fault maps are inherently incomplete; this might bias the performance criteria analysis. Furthermore, most faults form a topographic imprint, and therefore, integrating topography in the input data would be useful. Eventually, we will need to apply our new method to different datasets (different image types and resolutions, including satellite and aerial images), so as to verify that our novel approach can be generalized.

Figure 1 .

 1 Figure 1. Proposed conditional GAN architecture and methodology. The network has one Generator and one Discriminator sub-networks. The input of the Discriminator is the output of Generator combined with the input image. As training progresses, the Generator produces segmented images resembling more and more to the ground truth. The Discriminator then learns to discriminate between the ground truth image and the images produced by the Generator.

Figure 2 .

 2 Figure 2. Example of a convolution of an image X with a Gaussian blur kernel (or filter) K. X with size 8 × 8 passes through a convolution kernel function K with size 3 × 3. The filter moves toward the right of the image and new pixel values are computed. The sliding window (stride) is of 1 pixel.

Figure 3 .

 3 Figure 3. Optical image and manual expert annotations. (a) image of entire site; (b) annotations with major faults in red/pink and more minor faults in green; (c) patch of size 256 × 256 within input image; (d) manual annotations of the patch, with in red the major faults, in yellow the more minor faults, and in green the uncertain fault traces.

Figure 4 .

 4 Figure 4. An image patch, 1800 × 840, taken from the test region (a), the corresponding annotation made by human expert (major faults in red, more minor faults in yellow, and uncertain faults in black) (b) and our CGAN (c) prediction.

Table 1 .

 1 The confusion matrix to classify pixel values.

Table 2 shows

 2 IOU and Dice values for our proposed CGAN and the most classical U-net architecture, from the test images. The two performance indicators are better for the novel method we proposed than for U-net.

		IOU	Dice
	U-net	0.057	0.108
	Proposed CGAN	0.063	0.119
	Table 2. Performance of the proposed conditional GAN and
	U-net methods.