André Stumpf

Jean-Philippe Malet

Norman Kerle

Uwe Niethammer

Sabrina Rothmund

Image-based mapping of surface fissures for the investigation of landslide dynamics

Keywords: Landslide, Fissure, Geomorphological mapping, Remote sensing, Image filtering, UAV

The development of surface fissures is an important indicator for understanding and forecasting slope movements. Landslide investigations therefore frequently include the elaboration and interpretation of maps representing their spatial distribution, typically comprising intensive field work and instrumentation. Only recently aerial photography with sub-decimetre spatial resolution is becoming more commonly available and opens a window to analyse such features from a remote sensing perspective.

While these data are in principle helpful to elaborate maps from image interpretation techniques, there is still no image processing technique available to extract efficiently these geomorphological features. This work proposes a largely automated technique for the mapping of landslide surface fissures from very-high resolution aerial images. The processing chain includes the use of filtering algorithms and post-processing of the filtered images using object-oriented analysis. The accuracy of the resulting maps is assessed by comparisons with several expert maps in terms of affected area, fissure density and fissure orientation.

Under homogenous illumination conditions, true positive rates up to 65% and false positive rates generally below 10% are achieved. The resulting fissure maps provide sufficient detail to infer mechanical processes at the slope scale and to prioritize areas for more detailed ground investigations or monitoring.

Introduction

Observations of features and structural patterns of earth surface landforms can reveal information on the origin and mechanisms controlling the geomorphological processes. Structural geology and geomorphology have developed comprehensive concepts to delineate geomorphological units and structure types from remote sensing images, and infer about mechanical processes without necessarily measuring displacement, deformation or the applied forces directly [START_REF] Melton | Aerial photographs and structural geomorphology[END_REF][START_REF] Davis | Structural Geology of Rocks and Regions[END_REF][START_REF] Passchier | Microtectonics[END_REF][START_REF] Pollard | Fundamentals of Structural Geology, Third printing with corrections[END_REF]. Surface discontinuities observed in rocks and sediments have proven to be valuable indicators of the deformation history and stress pattern of the slope. For landslide analysis, their observation and interpretation can contribute to a better understanding of the controlling physical processes and help in the assessment of the related hazards [START_REF] Mccalpin | Preliminary age classification of landslides for inventory mapping[END_REF][START_REF] Fleming | Structures associated with strike-slip faults that bound landslide elements[END_REF][START_REF] Parise | Observation of surface features on an active landslide, and implications for understanding its history of movement[END_REF]. In hard-rock slopes, the analysis of structural discontinuities (faults, bedding planes, joints, fractures) allows to characterize potentially unstable areas [START_REF] Hoek | Rock Slope Engineering[END_REF][START_REF] Matheson | Rock Stability Assessment in Preliminary Site Investigations -Graphical Methods[END_REF][START_REF] Priest | Discontinuity Analysis for Rock Engineering[END_REF][START_REF] Selby | Hillslope Materials and Processes, 2nd Edition[END_REF][START_REF] Günther | Automated sliding susceptibility mapping of rock slopes[END_REF][START_REF] Jaboyedoff | Toward preliminary hazard assessment using DEM topographic analysis and simple mechanical modeling by means of sloping local base level[END_REF][START_REF] Glenn | Analysis of LiDARderived topographic information for characterizing and differentiating landslide morphology and activity[END_REF]. In soft-rock slopes and sediments, the analysis of surface fissures may indicate the development of future failures [START_REF] Krauskopf | Structural features of a landslide near Gilroy, California[END_REF][START_REF] Shreve | Sherman landslide, Alaska[END_REF][START_REF] Chowdhury | Tension cracks and slope failure[END_REF][START_REF] Abramson | Slope stability and stabilization methods[END_REF][START_REF] Khattak | Evolution of earthquake-triggered landslides in the Kashmir Himalaya, northern Pakistan[END_REF] and is often considered as a geo-indicator of the activity stage of a landslide. In sediments, the surface fissure characteristics also influence water infiltration and drainage, which in turn affect the ground-water system and the kinematic response of slopes to hydrological events [START_REF] Malet | Soil surface characteristics influence on infiltration in black marls: application to the Super-Sauze earth flow (southern Alps, France)[END_REF]Malet et al., 2005b;[START_REF] Van Asch | The diversity in hydrological triggering systems of landslides, The First Italian Workshop on Landslides[END_REF].

Maps of surface deformation features can be obtained by extensive field surveys either through the direct visual observation of the topography [START_REF] Fleming | Map and description of the active part of the slumgullion landslide[END_REF][START_REF] Meisina | Characterisation of weathered clayey soils responsible for shallow landslides[END_REF] or through the indirect measure of seismic wave propagation in tomography setups [START_REF] Grandjean | Characterisation of a landslide fissure pattern by integrating seismic azimuth tomography and geotechnical testing[END_REF][START_REF] Bièvre | Application of geophysical measurements for assessing the role of fissures in water infiltration within a clay landslide (Trièves area, French Alps)[END_REF].

Relatively large fissures on landslides may also be discernible in Very-High-Resolution (VHR) spaceborne images [START_REF] Glenn | Analysis of LiDARderived topographic information for characterizing and differentiating landslide morphology and activity[END_REF][START_REF] Youssef | Remote sensing applications to geological problems in Egypt: case study, slope instability investigation, Sharm El-Sheikh/Ras-Nasrani Area, Southern Sinai[END_REF] but typically those structures reach widths in the decimetre-range and at present only airborne photographs provide sufficient detail for their detection in the centimetric range. Recent studies [START_REF] Eisenbeiss | UAV Photogrammetry[END_REF]Niethammer et al., 2011a) have shown that VHR images acquired from Unmanned Aerial Vehicles (UAVs) are cost-efficient data sources for the monitoring of landslide surfaces with sub-decimetric image resolution. Especially small UAVs with payloads below 5 kg and operating altitudes below 2000 m are expected to be employed much more frequently in coming years [START_REF] Frost | Study Analysing the Current Activities in the Field of UAV[END_REF] though more specific regulations for their specific operational use are being discussed at national and international levels [START_REF] Prats | Regulations and Requirements[END_REF][START_REF] Watts | Unmanned Aircraft Systems in Remote Sensing and Scientific Research: Classification and Considerations of Use[END_REF].

Visual interpretation of VHR imagery is a classical method in geomorphology, but it remains subjective, and rather impractical for repetitive observations or the inspection of large areas. An increasing number of studies therefore targeted the development of automated techniques to extract relevant features from imagery [START_REF] Graham | Maximizing the accuracy of image-based surface sediment sampling techniques[END_REF][START_REF] Martha | Characterising spectral, spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods[END_REF][START_REF] Stumpf | Object-oriented mapping of landslides using Random Forests[END_REF]. Although, the detection and extraction of linear features is a fundamental operation in digital image processing [START_REF] Quackenbush | A review of techniques for extracting linear features from imagery[END_REF][START_REF] Mendonca | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF][START_REF] Papari | Edge and line oriented contour detection: State of the art[END_REF], relatively few studies have explored the application of automatic approaches for the mapping of geomorphological relevant linear features [START_REF] Baruch | Detection of gullies in roughly textured terrain using airborne laser scanning data[END_REF][START_REF] Shruthi | Object-based gully feature extraction using high spatial resolution imagery[END_REF].

Considering the increasingly widespread availability of sub-decimetre resolution images from UAVs and other airborne platforms, this study targeted the development of a semi-automatic image analysis technique to support geomorphologists in the detection, mapping and characterization of landslide surface fissures from VHR aerial images. In this context the term "semi-automatic" expresses that the technique requires user input to be adapted for different image types and environmental settings. The developed method is based on a combination of Gaussian directional filters, mathematical morphology and objectoriented image analysis (OOA) and was tested on a set of multi-temporal VHR images acquired at the Super-Sauze landslide (Southeast French Alps). The obtained results were compared to manual mappings carried out by experts combining image interpretation and field surveys.

Types of surface fissure observed on landslides

Detailed observations of landslide surface fissures were provided by [START_REF] Krauskopf | Structural features of a landslide near Gilroy, California[END_REF] who adapted analogies from structural geology for their interpretation and distinguished between strike-slip structures, normal faults, graben structures and compression structures. Also [START_REF] Ter-Stephanian | On the landslide cracks classification[END_REF] noticed the mechanical significance of surface fissures and elaborated a classification scheme relating fissure morphology and location within the landslide mass to corresponding mechanical processes. This included a first-order differentiation between upper extension, side friction, central compression, and lower creepon cracks. Although some authors used similar classification schemes [START_REF] Bombard | Une approche des problèmes posés par l'étude des mouvements de terrains[END_REF], the adopted terminology varies among different authors and affected lithologies [START_REF] Fleming | Structures associated with strike-slip faults that bound landslide elements[END_REF][START_REF] Cruden | Landslides Types and Processes[END_REF][START_REF] Fleming | Map and description of the active part of the slumgullion landslide[END_REF][START_REF] Walter | Joint analysis of the Super-Sauze (French Alps) mudslide by nanoseismic monitoring and UAV-based remote sensing[END_REF] and the terms crack and fissure are often used synonymously to refer to a variety of surface discontinuities.

Here, fissure is adopted as a generic term for open fractures on the topographic surface of a natural slope.

At first instance, transversal, longitudinal and diagonal fissures are distinguished according to their main orientation axes relative to the dip of the slope. This terminology can be used ad hoc to classify fissures solely based on geometric properties observed in the field or in an image. A more refined mechanical classification such as provided in [START_REF] Ter-Stephanian | On the landslide cracks classification[END_REF] will generally require considerations of the fissure patterns, the involved material and the local geometry of the slip surface. The term crack is used in this manuscript when referring to genetic processes described within classical fracture mechanics [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]. It should be noted that the term crack is also often adopted to refer to shrinkingswelling induced fractures [START_REF] Malet | Soil surface characteristics influence on infiltration in black marls: application to the Super-Sauze earth flow (southern Alps, France)[END_REF] which are not the objective of this study.

Classical fracture mechanics postulates tensile opening, sliding and tearing as the three basic modes for crack propagation (Fig. 1a). The concept has been developed for brittle material but is also adopted to explain fracturing of plastic materials at high deformation rates [START_REF] Schulson | Creep and Fracture of Ice[END_REF]. Surface fissures may develop from a combination of all three modes, whereas in practice, considering the relatively low tensile fracture toughness of most geomaterials [START_REF] Backers | Fracture Toughness Determination and Micromechanics of Rock Under Mode I and Mode II Loading[END_REF][START_REF] Ke | Determination of fracture toughness of anisotropic rocks by boundary element method[END_REF][START_REF] Schulson | Creep and Fracture of Ice[END_REF], tensile fracturing can be expected to dominate the formation of fissures at the free surfaces of a landslide. However, interpreting tension cracks as a direct indicator for a purely tensile stress regime may often fall too short. In fact, tensile fracturing may also result from relaxation of tensile stresses that originate from deformation induced by shearing and compression as well [START_REF] Wang | Brittle fracture in compression: Mechanisms, models and criteria[END_REF]. A mechanical interpretation and classification of the fissures must therefore consider the fissure pattern, the material and the landslide geometry.

FIG. 1 SOMEWHERE HERE

Fig. 1b-d illustrates three typical fissure patterns that are frequently used as geoindicators of specific deformation processes in the above-cited studies. One commonly observed example for such patterns is the formation of en-echelon fissure arrays (Fig. 1c), often also termed Riedel shears [START_REF] Riedel | Zur Mechanik geologischer Brucherscheinungen[END_REF]. They accommodate tensile stress and shear stress typically resulting from shear in the bounding zone of blocks moving with different displacement rates. Certain patterns such as arrays of transversal fissures (Fig. 1b) are typically associated with tension in the steeper upper slopes, whereas fissures resulting from compression and lateral extension (Fig. 1d) are more typically associated with gentler slopes in the transit and accumulation zones of landsides [START_REF] Sowers | Field Investigation[END_REF]. For landslides with a complex geometry, the position of those fissure patterns may however deviate considerably from this simple scheme (Niethammer et al., 2011a).

The Super-Sauze mudslide is an active slow-moving landslide located in the Barcelonnette Basin in the Southern French Alps (Fig. 2) that developed in weathered black marls in the 1960s, and features highly variable displacement rates (from 0.01 m.day -1 to 0.40 m.day -1 ) controlled by the local hydrological conditions (Malet et al., 2005b). The landslide measures 950 m from the main scarp to the toe, and is up to 150 m wide. The moving mass has a clay-rich matrix containing up to 30% coarse gravel as well as larger boulders and blocks (Malet et al., 2005b). The surface displays the signs of deformation in the form of ridges, bulges, lobes and fissures but also markers of surface erosion such as rills and small gullies.

Unlike surrounding stable areas the landslide surface is largely bare and only at a few locations, especially at its toe, cushion plants form small vegetation patches. Fissure widths of 0.01-0.40 m, lengths of more than 1.00 m and depths of up to 1.50 m [START_REF] Espinosa | Analysis and quantification of preferential flow on the Super-Sauze landslide[END_REF] can be observed in the field (Fig. 4b).

During the last 15 years, the landslide has been investigated through numerous monitoring campaigns including in-situ geophysical measurements, terrestrial and airborne LiDAR (Light-Detection and Ranging) and the acquisition of VHR optical imagery. In the VHR airborne optical images, the fissures can be recognized as dark curvilinear structures (Fig. 2c-e) as soon as their width approaches one pixel in size. Previous studies [START_REF] Malet | Les 'glissements de type écoulement' dans les marnes noires des Alpes du Sud. Morphologie, fonctionnement et modélisation hydro-mécanique[END_REF]Niethammer et al., 2011a;[START_REF] Walter | Seismic monitoring of rockfalls, slide quakes, and fissure development at the Super-Sauze mudslide, French Alps[END_REF] already discussed relationships between the observed fissure patterns (Fig. 2c-e) and strain resulting from a spatially heterogeneous displacement field and interactions between moving mass and the stable bedrock.

However, a full reconstruction of the complex bedrock geometry that may allow for a more detailed characterisation of the underlying deformation mechanisms has been completed only recently [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

FIG. 2 SOMEWHERE HERE

Airborne acquisitions of VHR optical imagery at the Super Sauze landslide

Between April 2007 and October 2009, diverse imaging systems and airborne platforms were used to acquire VHR images of the landslide at five different dates (Fig. 3). In July 2008, October 2008, and October 2009, a low-cost UAV system equipped with compact camera was operated at flight heights between 100 m and 250 m yielding images of the surface with a ground resolution between 0.03 m and 0.10 m. The individual images were corrected for barrel lens distortion, rectified according to ground control points (GCPs) measured with differential GPS (DGPS), and finally merged into one large orthomosaic. Further details on the image acquisition and processing were provided by [START_REF] Niethammer | UAV-based remote sensing of landslides[END_REF]2011b) who quantified the residual positional error (x-y) for the October 2008 images with 0.5 ± 0.57 m within the boundaries of the sliding area. The UAV images for July 2008 and October 2009 are expected to provide better positional accuracies because they were orthorectified using elevation models that were generated from a photogrammetric analysis of the images.

During the airborne LiDAR surveys (see § 3.2) of, respectively May 2007 and July 2009, two orthomosaics of optical images with full coverage of the landslide were recorded using medium format cameras (Fig. 3) mounted on, respectively an airplane and a helicopter. The surveys used fully integrated systems for direct georeferencing and orthorectification with LiDAR surface models (see § 3.2), which in general provide sub-decimetre positional accuracy in the x-y plane [START_REF] Vallet | GPS/IMU and LiDAR integration to aerial photogrammetry: Development and practical experiences with Helimap System®[END_REF].

For the study presented here, additionally 60 homologous tie points on stable areas were manually selected in the available images and showed a mean relative alignment error of 0.76 ± 0.82 m among the different acquisitions.

Further details on the adopted camera systems and the resolutions of the images resulting from the five surveys are summarized in Fig. 3. The figure also illustrates the considerable radiometric differences among the five images originating from illumination changes, seasonal variations and the distinct characteristics of the sensors. The scenes for May 2007, October 2008 and October 2009 were acquired under cloudy conditions with diffuse sky radiation and consequently show a more homogenous illumination of the surface. The scenes for July 2008 and July 2009 in contrast were recorded under sunny sky yielding strong contrast and many cast shadows. The latter are more prominent in the image for July 2008 which was recorded in the morning hours, at a relatively low sun angle. Available methods for absolute and relative radiometric correction can be employed for the radiometric alignment of satellite images [START_REF] Hong | A comparative study on radiometric normalization using high resolution satellite images[END_REF][START_REF] Vicente-Serrano | Assessment of radiometric correction techniques in analyzing vegetation variability and change using time series of Landsat images[END_REF] but, to the best of our knowledge, no approach exists to accurately align the radiometry of sub-decimetre images from different sensors, with substantial changes in illumination, a complex topography and changing surface characteristics. Initial test using histogram-matching, linear-regression [START_REF] Schott | Radiometric scene normalization using pseudoinvariant features[END_REF] and iteratively re-weighted regression [START_REF] Canty | Automatic radiometric normalization of multitemporal satellite imagery with the iteratively re-weighted MAD transformation[END_REF] did not provide satisfactory results. Consequently no radiometric normalization was performed and the image analysis technique was designed and tested with radiometric diverse imagery.

In order to calibrate adjustable parameters of the detection algorithm to the targeted fissures and the variable scene characteristics, the processing was first tested on a subset of the terrain covering ~ 14.000 m 2 in the central part of the landslide (Fig. 2a, Fig. 3). This section was characterized by different fissure patterns and recorded during all surveys (including July 2008 and October 2009 which did not yield full coverage of the surface). Subsequently, the developed workflow was applied on the full scenes for a comprehensive mapping and analysis of the fissure distribution. Corresponding results for the full extent of the Super-Sauze landslide and their mechanical significance are discussed in section § 5.2.

FIG. 3 SOMEWHERE HERE

LiDAR DTM

Two airborne LiDAR surveys were conducted in May 2007 and July 2009, respectively. The first survey used a Riegl LMS-Q560 laser scanner mounted on an airplane flying 600 m above the ground and resulted in a mean point density of 0.9 pts.m -2 after vegetation filtering. The residual 3D positional error of the ground points was quantified as 0.12 m. The second survey was conducted with a Riegl Q240i laser scanner mounted on a helicopter and after vegetation filtering resulted in a mean point cloud density of 3.2 pts.m -2 . The residual 3D positional error of the ground points was 0.07 m. Continuous surface raster with a pixel size of 0.5 m were interpolated from the respective point clouds using Delaunay triangulation. The resulting surface was then adopted for the extraction of the principal hydrological drainage lines.

Reference datasets: expert maps of surface fissures

Reference mappings of the fissure characteristics (type, distribution) were elaborated by an expert geomorphologist familiar with the study site. The fissures were first identified on-site during a field survey carried out in October 2009 at the same time as the acquisition of the UAV images. The position of the fissures was mapped using a dGPS survey and terrestrial photographs. Then image interpretation rules were defined to identify and digitize the fissures on the images as polyline vectors using a 2D view and at a scale of 1:250. The image interpretation rules were then applied to the four other images in order to elaborate an expert fissure map for each date. The resulting five maps were adopted as a reference to assess the performance of the semi-automatic method.

Image processing methods

While first generic edge detection operators were already proposed in the 1980s [START_REF] Marr | Theory of edge detection[END_REF][START_REF] Canny | A Computational Approach to Edge Detection[END_REF], the extraction of linear features from imagery remains a challenging task in many disciplines such as medical research [START_REF] Mendonca | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF], earth science [START_REF] Shao | Application of a fast linear feature detector to road extraction from remotely sensed imagery[END_REF][START_REF] Shruthi | Object-based gully feature extraction using high spatial resolution imagery[END_REF] or signal processing [START_REF] Lampert | A detailed investigation into low-level feature detection in spectrogram images[END_REF]. For our focus, the specific challenges posed for an automation of fissure detection can be summarized as follows:

 The approach should be scalable to apply for variable fissure sizes and image resolutions, and as insensitive as possible to variable radiometric image characteristics;  The technique should not respond to edges but enable the detection of dark curvilinear structures that may be oriented at any direction. Classical techniques such as Sobel operator and the Canny detectors [START_REF] González | Digital Image Processing[END_REF] have been designed specifically for edge detection and are not directly applicable;

 The complex micro-topography, the presence of rock blocks and gravels as well as small patches of vegetation yield highly textured images. Consequently, the approach should enable to smooth out spurious signals from the noisy background while still retaining small partially disconnected linear features of interest. Contextual scene information should be taken into account to resolve ambiguities of the local features.

Considering these challenges, a processing workflow including three main stages was developed. Firstly, a set of scalable Gaussian filters is applied to detect fissure candidates and suppress responses at edges. Secondly, a set of morphological filters is used to close small gaps along the extracted candidates.

Thirdly, an object-oriented procedure is followed to eliminate some of the false positives exploiting higher-level scene information with contextual rules.

FIG. 4 SOMEWHERE HERE

Stage 1: Extraction of fissure candidates using a Gaussian matched filtering algorithm

A particularly well-studied example for the detection of dark curvilinear structures is the extraction of dark blood vessels in photographs of the human retina. Based on the observation that the cross-profiles of the vessels resembles a Gaussian distribution, [START_REF] Chaudhuri | Detection of blood vessels in retinal images using two-dimensional matched filters[END_REF] proposed the use of a matched filter (MF) that is essentially a Gaussian convolution kernel subtracted by its own mean value. As illustrated in Fig. 4a, the cross-sections of surface fissures can be approximated with a Gaussian distribution and an MF scaled to the size of the fissure will give a peak response when crossing the fissure at an angle of approximately 90°. Because the MF still yields errors such as false detections at step edges (Fig. 5a,c) numerous extensions [START_REF] Hoover | Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response[END_REF][START_REF] Sofka | Retinal vessel centerline extraction using multiscale matched filters, confidence and edge measures[END_REF] and alternative approaches [START_REF] Mendonca | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF][START_REF] Soares | Retinal vessel segmentation using the 2-D Gabor wavelet and supervised classification[END_REF] have been developed. Recently, [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF] proposed modification to the original MF filtering approach integrating a first order derivative of a Gaussian function (FDOG) to locally adapt the thresholds separating dark line from non-target features.

Compared to other state-of-the-art algorithms their approach provided competitive accuracies while being a computationally efficient and hence easier to apply on the large images resulting from VHR remote sensing.

For this study, a similar approach was implemented in ENVI-IDL 4.8 (ITT Visual Information Solutions).

The algorithm and its parameterization are detailed below.

FIG. 5 SOMEWHERE HERE

The MF is a two dimensional kernel defined in the x-direction by an inverted Gaussian profile (Fig. 5b), and in the y-direction by replicates of the same profile (Fig. 5d). It may be denoted as:

, ; √ , for |x|≤3σ, |y|≤L/2
Eq. 1

where σ denotes the standard deviation of the Gaussian functions and relates to the width of the targeted feature. To centre the kernel on zero, it is subtracted by its own mean m. The extent of the kernel in the xdirection is typically constrained to 3σ, whereas L defines the extent of the kernel in the y-direction and can be related to the length of the fissures. Because the matched filter still yields false responses at dark and bright step edges (Fig. 5c) [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF] proposed to use the response of the FDOG to locally adjust the thresholds which are applied to classify the MF response into fissure and non-fissure structures.

In analogy to Eq. 1, the first order derivative filter may be denoted as:

′ , ; √ , for |x| 3σ, |y| /2
Eq. 2 Fig. 5f illustrates that the FDOG responds with a single peak to edges but with a zero crossing at the centre of the idealized fissure. A simple mean filter can be applied to broaden the zero crossing into a plateau covering the whole width of the fissure (Fig. 5f). Subtracting the smoothed FDOG response from the MF response will attenuate the signal at edges while at the position of the fissure the full response is retained (Fig. 5h).

Since the orientation of the fissures is a priori unknown, multiple rotated versions of the Gaussian filters are applied on the image and for each pixel only the maximum response value is retained. This corresponds to finding the angle , which maximizes the filter response at a given position in the image , using Eq. 3:

, arg max , ⊗ , for 0 Eq. 3

where ⊗ denotes the convolution operator and the orientation of the MF.

The calculation of the maximum response image R can then be obtained with Eq. 4: where M denotes the above-mentioned mean filter used to broaden the zero crossing to the width of the fissures.

While [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF] used a very broad mean filter with a fixed size, we suggest to use a kernel size that matches the width of the Gaussian kernel (6σ) and is thereby related to the width of the targeted features (Fig. 5a,f). In contrast to early studies where the FDOG response was used to locally adapt the threshold [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF] the final response image is obtained by subtracting the FDOG from the GMF response using Eq. 6:

, , * , Eq. 6

where Ct denotes a user defined trade off parameter to adjust the sensitivity of the detection with typical range of values between 3 and 4. A threshold T is defined by Eq. 7:

2 Eq. 7

where is the mean of the response image R and .

A binary fissure candidate map F map is obtained by applying the threshold T on the response image R using Eq. 8:

, , ∶ 1 and

, , ∶ 0
Eq. 8

The thresholding after subtraction of the FDOG response was found to provide a generally more robust attenuation of undesired edge responses than the technique previously applied by [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF].

In summary, the user needs to specify four simple parameters, namely (1) the scale of the filter kernels in terms of σ, (2) the length L of the kernel, (3) the constant Ct of the thresholding sensitivity and ( 4) the number of orientations n θ at which the filters are calculated. In this study, n θ was kept constant at 36 for all experiments, whereas, if computational time becomes an issue, the angular resolution may be reduced to 12 steps without major losses of accuracy. To determine σ a tool was created,which allows drawing profiles on the image and automatically estimates the fitting Gaussian function (Fig. 4). Cross-profiles of the smallest fissures visible in the image with the coarsest resolution (0.1 m.pixel -1 resolution) were best fitted by Gaussian curves with σ ≈ 0.6. To ensure a homogenous scale of the detected features among all images, the kernel can be scaled by changing σ relative to the image resolution. If, for instance, the image resolution is increased to 0.08 m.pixel -1 , a value of σ ≈ 0.75 yields a kernel with the same physical size (Table 1). The same applies for the filter length L which was estimated as 1 m corresponding to the typical minimum length of the fissures. Resampling of the images can thereby be avoided. In our experience, σ establishes the lower bound for the width of the targeted features, whereas the filters still remain sensitive to features which are up to 5 times larger. For the choice of σ it is also helpful to note that the discrete kernel cannot represent FDOG functions with σ ≤ 0.5.

To assess the sensitivity of the parameters and to determine a suitable threshold parameter Ct, a sensitivity analysis was carried out on a subset of the October 2008 image. Based on a visual assessment, values of L = 1 m and σ = 0.75 were found suitable for the detection of the fine fissure structures. The preliminary analysis also showed that increasing the parameters L and σ directs the detection towards more elongated and broader features, whereas in general the sensitivity of those parameters is rather low compared to the influence of the threshold Ct. Values of Ct = {0.0, 1.0, 2.0, 3.0, 4.0} were tested and based on a visual assessment of the outputs a value of Ct = 3 was established for an optimal trade-off between detection rate and the amount of false positives. The final parameter set is summarized in Table 1.

TABLE 1 SOMEWHERE HERE

Stage 2: Connection of broken lines using structuring elements

The highly textured surface of the landslide constitutes a noisy background that affects the detection especially at section where the fissures are very thin or partially occluded. While a human operator can easily interpolate broken lines through perceptual grouping [START_REF] Metzger | Gesetze des Sehens[END_REF], this needs special attention for a semi-automated mapping technique.

To close small gaps between broken line segments of the detected candidates, a hit-or-miss transform algorithm [START_REF] Serra | Image Analysis and Mathematical Morpholgy[END_REF] was used. The transform assigns a value of 1 to each pixel whose local neighbourhood fulfils the criteria defined by hit-and miss structures (Fig. 6a), also known as structuring elements. They were defined to address all plausible 3-by-3 neighbourhoods representing small gaps in the detection starting from four prototype hit-structures shown in Fig. 6b. The respective miss-structures (Fig. 6c) are typically derived by simply inverting the prototype hit-structures, and both elements were rotated (Fig. 6d) to test for a total number of 24 possible neighbourhood arrangements. Exceptional cases were thereby the structuring elements for closing directly diagonal gaps, where an extended neighbourhood was used for the hit-and miss structures (Fig. 6b,c) to prevent connections of lines running parallel to each other.

FIG. 6 SOMEWHERE HERE

The connectivity of the line segments was also particularly important for the subsequent object-oriented post-processing, where objects constitute from pixels groups connected in a Von Neumann neighbourhood (4 adjacent pixels at each side), and small isolated objects could be disregarded as noise.

Stage 3: Object-oriented analysis for false positive removal

Due to visually similar objects, such as linear erosion features (rills, small gullies) and elongated shadows induced by the micro-topography, the fissure candidates resulting from the described filtering routine may still comprise numerous false positive detections. While a human interpreter can differentiate most of the false positives assessing the geospatial context of the scene, the efficient use of such information with automated systems is a challenge for object-oriented image analysis. To exploit the contextual scene information for an automatized refinement of the extracted fissure candidate maps, an object-oriented routine that integrates spatial reasoning into an explicit form was elaborated and implemented using eCognition 8.64 [START_REF] Trimble | eCognition Developer 8.64.0 Reference Book[END_REF]. The elaborated routine included the following steps:

1) The ratio of shadow around the detections is evaluated and candidates with a ratio of shadow pixels in their smallest enclosing circle above 33% are regarded as false detections induced by shadings of the micro-topography. This ratio threshold was determined empirically through visual inspection of the candidate fissures, and selected to capture elongated false detections with one side lying fully in shaded zones. The threshold for shadow can thereby be adjusted according to the illumination conditions and the dynamic range of the image (Table 2).

2) Further false detections may result from vegetation which typically shows a lower reflectance in the green and red channel compared to the blue. The blue ratio in the sum of all channels is consequently typically below one third for vegetated areas. The suitable value varies slightly with the illumination conditions and the season, and Otsu's method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF] was employed to automatically adapt to such changes. Through an iterative testing of all possible values, Otsu's method determines threshold value that maximizes variance between two classes in an image. Hence, constraining the search space to all pixels with a ratio blue below 33%, the algorithm was used to determine the thresholds that maximize the contrast between vegetation and the background (Fig. 7). Fissure candidates covered by the resulting vegetation class, or having a relative border length larger than 0.15, were subsequently removed.

FIG. 7 SOMEWHERE HERE

3) Another class of frequent false detections resulted from linear objects such as rills, gullies and nearly vertical steps at the landslide flanks, which may locally obtain similar characteristics as the targeted fissures. To test for the presence of larger linear features and evaluate their relationship with fissure candidates, a strategy to suppress additional false positives was required. For the mapping of the larger linear elements two sources were adopted. First drainage lines were extracted from the LiDAR DTMs using hydrological standard tools [START_REF] Tarboton | On the extraction of channel networks from digital elevation data[END_REF] and enlarged with a surrounding buffer of 0.5 m. A second approach was to repeat the Gaussian filtering with the parameter set indicated in Table 1, but with a two times increased scale σ and a 5 times coarser image resolution (resampled with bilinear interpolation). This is equivalent to a search with a 10 times larger scale providing a sufficiently large scale difference to assure that none of the detected linear features would correspond to fissures. The linear objects extracted with both approaches were virtually overlaid with the fissure candidates, and the difference of the orientations of their respective centre lines was adopted as criteria to evaluate if the fissure candidate was in fact part of a larger linear object or constitutes an independent structure (Fig. 8). Image-based measurements of the angular offset of the fissured indicated a minimum offset of about ± 13°. Considering that the lowest effective friction angle values measured for the landslide material are ' = 26° (Malet et al., 2005a), the thresholds are consistent with the orientation of '/2 that the Coulomb criterion predicts for the orientation of shear fissures at the landslide boundary [START_REF] Tchalenko | Similarities between Shear Zones of Different Magnitudes[END_REF].

FIG. 8 SOMEWHERE HERE

4) A last filtering step was implemented by removing all candidates with length not longer than 0.4 m and an area smaller than 0.1 m 2 . Finally all fissure candidates falling in areas with a fissure class density lower than 1% in a surrounding neighbourhood of 10 m 2 were considered as noise and also removed.

Table 2 displays that most adopted thresholds were kept the same among all the images and only the classification rule for the shadow areas was adapted in order to compensate radiometric differences in the input images.

TABLE 2 SOMEWHERE HERE

Results

Comparison with multi-temporal manual mappings

The primary output of the developed processing routine is a map of the detected fissures represented by polygons. Applying a Delaunay triangulation that extracts the skeleton of those polygons [START_REF] Trimble | eCognition Developer 8.64.0 Reference Book[END_REF], a 2D line representation, which enables a more immediate comparison with expert mappings, can be obtained.

FIG. 9 SOMEWHERE HERE

Fig. 9 displays an example of comparison between an expert map and the result of the semi-automatic detection. A first visual assessment of the obtained maps suggested better agreement of the fissure patterns in areas with high contrast and low texture (Fig. 9a), whereas false positives and false negatives concentrated in sections with low contrast and increased surface texture (Fig. 9b).

For a quantitative assessment of the mapping accuracy the obtained results were compared with the expert mappings in the central part of the landslide (Fig. 9c) at all five dates. While several accuracy measures for geographic line datasets have been already proposed, there is still no consensus about one generally applicable technique and the metrics should be selected according to the problem at hand (Ariza-López and Mozas-Calvache, 2012). Here, we focused on three crucial aspects of the map accuracy that may have direct implications for their further use, namely the size of the affected (e.g. fissured) area, the length and density of the fissures, and their orientation.

Size of the fissured area

Tveite and Langaas (1999) suggested an accuracy measure for line datasets based on repeated buffering and overlay operations of detected and reference line datasets. A similar strategy was adopted in this study by repeatedly calculating true positive and false positives rates from two raster representing the detections and the expert mapping at increasingly coarser resolutions. Raster were calculated at 10 cm steps for resolutions between 0.10 m and 1.00 m, and each pixel was assigned as fissured or non-fissured area according to the presence or absence of a fissure in the detections and the reference map, respectively. The resulting Receiver Operating Characteristics (ROC) plots are presented (Fig. 10). The analysis showed a correspondence with the expert maps at true positive rates typically above 40% and up to 65%. The false positive rates were below 5% except for the scenes recorded with full sunlight where false positive rates up to 9% could be observed (Fig. 10).

FIG. 10 SOMEWHERE HERE

Fissure length and density

Hydrological models that integrate the influence of surface fissures on infiltration and preferential flow have demonstrated that the fraction of fissures per unit area is an important parameter with considerable influence on the modelled water storage (Malet et al., 2005b;[START_REF] Krzeminska | A conceptual model of the hydrological influence of fissures on landslide activity[END_REF]. Such models are typically generated at slope scale with grid resolutions below 10 m. To assess the accuracy of the extracted maps with respect to this potential application, the fissure density was calculated as the line length in circular sliding windows with diameters between 2 and 10 m, and compared among automated detection and expert mappings.

FIG. 11 SOMEWHERE HERE

The regression plots in Fig. 11 illustrate the correlation of the fissure density estimates with a 5 m circular sliding window yielding coefficient of determination (R 2 ) typically above 0.5. The regression analysis further indicated generally higher densities resulting from the semi-automatic detection originating from false positive detections but also from a stronger generalization of the fissure line drawings within the expert mapping. Exceptions from this general trend are the results obtained from the image of July 2008

which was recorded at a low sun incidence angle leading to a relatively low = 0.36. The bar plots in Fig. 11display the generally higher R 2 values at increasing resolutions of the density raster. This is a wellknown effect of spatial aggregation on correlation statistics [START_REF] Gotway | Combining incompatible spatial data[END_REF] but also reflects the contrast between stronger discrepancies of local details and a better correspondence of the global fissure pattern pictured in the respective maps. The highest correlation was observed among the mappings for May 2007 with 0.88 indicating that the lower resolution of the corresponding input image was not an important factor for the accuracy of the detection.

Fissure orientation

As outlined in the introductory section, different fissure patterns may signal respective mechanical processes and statistics of the principal fracture orientation often allow to estimate the directions of the principal stresses [START_REF] Pollard | Fundamentals of Structural Geology, Third printing with corrections[END_REF]. The fissure orientations were quantified as a third factor to assess the accuracy of the extracted maps using rose diagrams frequently employed for the analysis and interpretation of two dimensional orientation data [START_REF] Jammalamadaka | Topics in Circular Statistics[END_REF]. Rose diagrams with a bin width of 10° were computed on a 10 m regular grid for the semi-automatic detections and the expert mappings at all five dates. Considering the length and direction of each bin expressed as a respective vector the preferred fissure orientation within a grid cell can be calculated by summing the vectors over all bins. Taking into account all cells containing fissures in both the expert map and the semi-automatic detection, the Mean Absolute Error (MAE) of the mean orientations provides a quantitative measure for the orientation accuracy.

The rose diagrams plots and error statistics in Fig. 12 with displacement rates significantly above the average annual rates [START_REF] Travelletti | Imagerie multi-paramètres et multi-résolutions pour l'observation et la caractérisation des mécanismes de glissements-coulées[END_REF].

Comparing the detection results of May 2007 (Fig. 13a) and October 2008 (Fig. 13b), a significant increase in the abundance of fissures could be noted for the entire landslide. This can be attributed to a phase of strongly increased displacement rates (up to 3.5 m.day -1 ) in early June 2008 [START_REF] Travelletti | Imagerie multi-paramètres et multi-résolutions pour l'observation et la caractérisation des mécanismes de glissements-coulées[END_REF] preceding the UAV survey in October 2008. Though, in October 2008, the displacement rates already consolidated again at average rates between 0.01 and 0.03 m.day -1 , most of the fissures induced in June were preserved and evolved at the surface until October. This view is supported by the results obtained for the test area with the five scenes (Fig. 10, Fig. 11) picturing rather a transient evolution than a complete reorganization of the fissure patterns. Despite partially strong disagreement in the absolute measured fissure density, both expert maps and semi-automatic mapping showed an increase in fissure density after May 2007, with higher values in October 2008 (Fig. 11) than directly after the peak displacement in spring. Pluviometric records for the area in 2008 show relatively dry summer season with cumulative rainfall of 110 mm for the month of July, August and September, suggesting that the increased fissure density in October is partially caused by an increased brittleness of the upper soil layer that dried out during summer.

Besides the general increase in the amount of fissures, it is intriguing to observe that at several local plots, similar fissure patterns can be observed at approximately the same positions through time (Fig. 12, Fig. 13a,b), despite maximal displacements of up to 55 m between October 2008 and October 2009 (Niethammer et al., 2011a). This indicates the recurrent continuous in-situ formation where the fissures provide a close representation of the local strain field, similar as observed for the evolution of glacier crevasses [START_REF] Harper | Crevasse patterns and the strain-rate tensor: a highresolution comparison[END_REF].

Previous studies [START_REF] Malet | Les 'glissements de type écoulement' dans les marnes noires des Alpes du Sud. Morphologie, fonctionnement et modélisation hydro-mécanique[END_REF]Niethammer et al., 2011a;[START_REF] Walter | Seismic monitoring of rockfalls, slide quakes, and fissure development at the Super-Sauze mudslide, French Alps[END_REF] already observed close relationships between the occurrence of fissures and the geometry of the stable bedrock at the Super-Sauze landslide. They also noted a general contrast between higher water content and rather ductile behaviour in the lower subsurface (< 1m) and typically lower water content of the topsoil yielding more brittle behaviour at the surface. The surface fissures can therefore be understood as the response to stresses induced in the topsoil through coupling with ductile strain in the deeper subsurface. A similar model was already described by [START_REF] Fleming | Structures associated with strike-slip faults that bound landslide elements[END_REF] and adopted as a basis to qualitatively estimate the patterns of flow and stresses from a joint-interpretation (Fig. 13) of the detected fissure patterns and a geometrical model of the stable bedrock [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

Considering the bedrock geometry and a formation of the open fissures normal to the direction of the least compressive stress [START_REF] Pollard | Fundamentals of Structural Geology, Third printing with corrections[END_REF], three different flow field patterns leading to the fissure formation at the Super-Sauze landslide can be suggested. First, lateral shear at external and internal landslide boundaries aligned with the general flow field leads to the formation of diagonal shear fissure arrays (Fig. 13e). Second, longitudinal compressive and tensile strain related to abrupt changes in the slope of the sliding surface induces tensile stresses at the surface that results in transversal fissure arrays (Fig. 13c). Third, divergence of the flow field over topographic ridges and at the outlets of confining topographic channels induces lateral extension and tensile stresses resulting in longitudinal fissure arrays (Fig. 13d). At several locations, those processes overlap and lead to the formation of mixed structures such as a radial fissure patterns displayed in Fig. 4 and Fig. 13c, resulting from lateral shear and longitudinal strain, and from local divergent stress field, respectively.

FIG. 13 SOMEWHERE HERE 6 Discussion

Deformation patterns at the surface of landslides are important indicators for the mechanical processes, whereas the elaboration of detailed maps of such features remains a challenging and time-consuming task.

While [START_REF] Sowers | Field Investigation[END_REF] still argued that aerial photographs do not reach sufficient resolution for such mappings, modern digital sensors and new aerial platforms such as UAVs today provide the necessary level of detail. Furthermore, this study demonstrated the possible use of a semi-automatic image processing chain for the extraction of surface fissure maps.

The accuracy of the method was assessed by comparisons with expert maps and demonstrated heterogeneous areal accuracies with true positive rates of up to 65% and false positive rates generally below 10%. Also, the orientation accuracy showed a variable quality of the resulting maps with mean deviations between 9.7° and 22.5°. The fissure densities derived from both maps have significant correlations (R 2 between 0.36-0.78), whereas the semi-automatic detections yield typically higher estimates. Interestingly, this difference is more pronounced with the images of 2009 (Fig. 10) reflecting the contrast between increased semi-automatic detection rates at higher resolutions and the fixed scale of the expert mapping. Contrariwise, the best agreement among detection and expert maps was measured for the scene of May 2007 showing that the lower resolution does not necessarily yield lower accuracies.

Generally lower accuracies were observed for the scenes recorded with full sunlight at the surface in July In the initial stage of the processing chain, a low-level linear feature detector is used. Similar techniques yield competitive results in medical image analysis [START_REF] Zhang | Retinal vessel extraction by matched filter with firstorder derivative of Gaussian[END_REF], whereas the accuracies achieved with aerial images in this study are still significantly lower. This must be attributed to the generally higher complexity of outdoor scenes and at the moment still requires additional steps and parameters to take the contextual scene information into account. The use of an OOA heuristic-based post-processing technique proofed useful for the removal of false positives and helped to objectify the image analysis by transferring expert knowledge in an explicit form. The analysis still relies on a number of fixed thresholds which may hinder an easy transfer of the entire processing chain to a different geographic area. This concerns especially parameters that require knowledge of the local processes (e.g. minimum fissure length, effective friction angle) while thresholds that can be determined directly from the image (e.g. shadows, vegetation) may be adjusted more effortless.

The development of surface fissures precedes and accompanies especially slow-and very slow-moving landslides [START_REF] Cruden | Landslides Types and Processes[END_REF] making the developed technique particularly applicable to such types of landslides and to potentially unstable slopes. However, the spatial resolution of the acquired images must at least match, or should ideally exceed, the width of the targeted fissures, and the vegetation must be sufficiently sparse to permit direct view on the bare ground. The results of this study demonstrated that if those requirements are met, the obtained fissure maps can already provide sufficient accuracy to infer the landslide dynamics and mechanical processes (see § 5.2.) at the slope scale. Density maps from both semi-automatic and expert mappings show a strong spatial and temporal variability of the fissure abundance pointing toward important local and temporal contrasts in the infiltration capacity which may need considerations in the design of hydro-mechanical models. An analysis of the evolution and mechanics of individual fissures will however require higher temporal resolution and terrestrial cameras have recently been installed at the surface of the landslide to record imagery for further research in this direction. It would also be desirable to test the developed technique for the investigation of other landslides with different characteristics in order to validate a more general applicability of the approach and the mechanical interpretation of the observed fissure dynamics.

The OOA heuristics already considers multi-scale information to some degree (see § 4.33), whereas for further methodological improvements an explicit integration of an automatic scale selection technique at the low-level filtering stage appears as a promising approach to further reduce heuristics and tuneable parameters [START_REF] Stumpf | Multi-scale line detection for landslide fissure mapping[END_REF]. The first and second stages of the proposed method are generic for the detection of dark linear features and could in principle also be applied to detect other geomorphological and geological structures with similar characteristics. The proposed technique might be of interest for the mapping of erosion gullies [START_REF] Shruthi | Object-based gully feature extraction using high spatial resolution imagery[END_REF], geological lineaments [START_REF] Mallast | Semi-automatic extraction of lineaments from remote sensing data and the derivation of groundwater flow-paths[END_REF], ice-glacier crevasses [START_REF] Vaughan | Relating the occurrence of crevasses to surface strain rates[END_REF] or tectonically induced fractures [START_REF] Allmendinger | Invited review paper: Neogene to Quaternary tectonics of the coastal Cordillera, northern Chile[END_REF], sufficiently larger to be depicted in sub-metre satellite images.

Considering the intrinsic disagreement in expert mappings of linear features, especially in the inter-and extrapolation of lines [START_REF] Sander | Ground-water exploration based on lineament analysis and reproducibility tests[END_REF], further studies should also include an assessment of the uncertainties of reference maps since their quality can strongly bias the evaluation of different alternative approaches (Lampert et al., submitted).

Conclusions

This study developed an image processing chain to extract surface fissures from heterogeneous sets of VHR aerial images and tested the approach with a challenging multi-temporal set of images recorded at the Super-Sauze landslide for five different dates. The first two stages of the developed workflow combine families of Gaussian matched filters and morphological filters and are followed by an objectoriented analysis to reduce the amount of false positive detection using contextual information and auxiliary topographic information. The detection results can be represented in raster maps or optional by centre skeleton lines.

Under homogenous illumination conditions a comparison of the results with expert mapping demonstrated detection rates of up to 65% and orientation errors below 10°. Contrary, the technique is relatively sensitive to shading effects at full sunlight and prone to errors especially at low sun incidence angle.

A joint-interpretation of obtained fissure maps and of a 3D geometrical model of the stable bedrock demonstrated their complementary use for a better understanding of the geomorphological and geomechanical processes, such that the detected fissure pattern may be used for first approximation for mechanical processes in the recent deformation history of a slope.

Possible directions for further research are the reduction of tuneable parameters and a more immediate exploitation of multi-scale information, as well as an adaption of the technique to other linear features with geomorphological and geological relevance. [START_REF] Sowers | Field Investigation[END_REF]. [START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

Table 1

Parameter set of the Gaussian filters scaled according to the respective image resolution.

Table 2

Summary of the thresholds adopted in the object-oriented post-processing routine. See text for details.

  response values are automatically set to zero and only values greater than zero are retained. The FDOG filter is rotated according to the determined , and the corresponding response image D can be derived by Eq. 5:

  FIG. 12 SOMEWHERE HERE

  2008 and 2009 and the worst results were obtained for July 2008 when images were recorded at a relatively low sun incidence angle. Since the direct sunlight induces shading that affects the local contrast and global image normalization methods cannot alleviate this problem and image acquisition with diffuse skylight appears to be the generally better option.

Fig. 1 .

 1 Fig. 1. Generic types of surface fissures and their typical spatial occurrence within a landslide mass. (a) Modes of fracture propagation: mode I (opening), II (sliding) and III (tearing). (b) Fissures developing predominately in mode I and resulting from tensile stress. (c) Fissures developing predominately in mode I and resulting from shear stress. (d) Fissures developing predominantly in mode I resulting from compressive stress and lateral expansion. The central figure is adapted from[START_REF] Sowers | Field Investigation[END_REF].

Fig. 2 .

 2 Fig. 2. Oblique view of the Super-Sauze landslide combining a hillshade image derived from an airborne LiDAR DTM (July 2009) and a UAV image (October 2008). (a) The main scarp (hashed black line), transport and accumulation zone (black outline), and the area of interest for the multi-temporal analysis (white square) are indicated. UAV image subsets show (b) compression ridges, (c) longitudinal fissures, (d) diagonal fissures at the boundary of the active part, and (e) and transversal fissures.

Fig. 3 .

 3 Fig. 3. Subsets of orthophotographs (see location in the white bounding box in Fig. 2a) acquired at five different dates with details of the acquisition systems and image ground resolutions.

Fig. 4 .

 4 Fig. 4. (a) Subset (see extent in Fig. 3) of the UAV image from October 2008 showing typical fissure patterns and (I-IV) grey-value profiles (green channel) approximated with Gaussian curves. (b) Field terrestrial photograph taken in October 2009.

Fig. 5 .Fig. 6 .Fig. 7 .

 567 Fig. 5. Illustration of the principles of the Gaussian filtering for (a-c, f-h) a simplified 1-D case, (d, i) a 3D visualization of 2D filters and (e, j) the filter responses for the image subset in Fig. 4a. See text for details.

Fig. 8 .

 8 Fig. 8. Illustration of the object-oriented post-processing routine. (a) Fissure candidates that overlapped with linear structures. (b) Linear structure detected at a ten times greater filter scale. The fissure candidates aligned with the linear structures at angles below ± 13° were removed.

Fig. 9 .

 9 Fig. 9. Example of comparison of the obtained fissure maps with the expert mapping for October 2010 (fissures in red). (a) Area with relatively high agreement of the mapped fissure patterns. (b) Area with relatively high rate of false negatives and false positives. The scale of the representations corresponds approximately to the scale used for the expert mapping (1:250).

Fig. 10 :

 10 Fig. 10: Receiver operating characteristics (ROC) plots for the fissured area at different map resolutions.The sky conditions for the five different dates are indicated.

Fig. 11 .

 11 Fig. 11. Correlation between fissure density estimates at 5 m raster resolution based on semi-automatic detections and expert mappings from the five images. The bar plots at the bottom display the R 2 coefficient at different raster resolutions.

Fig. 12 .

 12 Fig. 12. Rose diagram plots with mean orientation (red line) and error statistics for the mean fissure orientation per 10 m grid cell for the test area at the five different dates. For visualization, the rose diagrams where plotted over a hillshade of the landslide surface and the scatterplot angles were centred at 90°.

Fig. 13 .

 13 Fig. 13. Pseudo 3D view showing the landslide dynamics inferred from the fissure patterns detected in the aerial images of (a) May 2007 and (b) October 2008. (c, d, e) Close up views for October 2008 showing inferred landslide dynamics and stress vectors. The results are overlaid on a hillshade model of the topography of the stable bedrock proposed by[START_REF] Travelletti | Characterization of the 3D geometry of flow-like landslides: A methodology based on the integration of heterogeneous multi-source data[END_REF].

Acknowledgements

The work was supported by the project SafeLand "Living with landslide risk in Europe: assessment, effects of global change, and risk management strategies" (Grant Agreement No. 226479) funded by the


Willy Allègre

Thomas Burger

Pascal Berruet

Jean-Yves Antoine

A Non-Intrusive Monitoring System for Ambient Assisted Living Service Delivery

Keywords: Pervasive Systems, Ambient Assisted Living, Model-Driven Engineering, Ontologies, Home Automation Services, Context-Awareness

Automation of smart home for ambient assisted living is currently based on a widespread use of sensors. In this paper, we propose a monitoring system based on the semantic analysis of home automation logs (user requests). Our goal is to replace as many sensors as possible by using advanced tools to infer information usually sensored. To take up this challenge, an ontology, automatically derived from a model-driven process, firstly defines user-system interactions. Then, the use of rules allows an inference engine to deduce user location and intention leading to adapted service delivery.

Introduction

Care for dependents is becoming a major social and economic issue for the next few years. In 2050, 30% of people from the european countries will be at least 65 years old [1], thus increasing the number of dependent people with chronic diseases. The number of potential caregivers can not evolve accordingly, so that home automation for ageing and disabled people is considered as an alternative solution to fill the gap, allowing one to facilitate or automate devices activation and to provide services tailored to the user needs: we talk about Smart Home.

Due to the miniaturization of electronic devices and their dissemination in the environment, smart home is no longer confined to the simple control of household appliances. This field of application is evolving more and more toward Pervasive Systems [START_REF] Weiser | Some computer science issues in ubiquitous computing[END_REF], where context-awareness is an important challenge. Indeed, knowing who is doing what and where is essential to adapt a "smart behavior", especially in the area of Ambient Assisted Living. The use of sensors distributed in the environment is the most commonly used solution to collect these contextual data. Two main methods can be distinguished, which can be coupled [START_REF] Yu | Human localization via multi-cameras and floor sensors in smart home[END_REF][START_REF] Nourizadeh | Medical and home automation sensor networks for senior citizens telehomecare[END_REF][START_REF] Long | Using graphs to improve activity prediction in smart environments based on motion sensor data[END_REF]: the use of cameras which is intrusive, and the massive deployment of sensors in the environment which implies a significant financial cost.

Previous works from [START_REF] Truong | Proactive remote healthcare based on multimedia and home automation services[END_REF] highlighted the interest of monitoring activities based on home automation logs (user requests only) to limit the use of these sensors. Their monitoring system is based on a frequency analysis in which a service is considered as a repeating pattern over time. Based on the assumption of the existence of very regular activities, the goal was to automatically extract and to propose daily living scenarios to wheelchair-bound users with limited capabilities. However, this assumption is not always verified in practice, so that some potential scenarios are not proposed. For example, the user needs to switch on/off lights in rooms leading to the bathroom at night. This kind of interaction can be requested from time to time, in a non-regular manner.

In this paper, instead of a frequency analysis, we propose to use a semantic analysis of each home automation log to interpret service activation at a given time. In this perspective, it is possible to provide the user adapted services without any training steps and using as few sensors as possible. To take up this challenge, we first propose to model user-system interactions and their semantics by means of an ontology (section 2). Then, the use of rules applied to ontology concepts and properties allows an inference engine to deduce user location and intention, leading to service delivery (section 3). Finally, in the context of a large scale deployment, we propose to derive automatically the ontology from a model-driven process developped in previous works (section 4). A real case study (section 5) presents the benefits of the entire work.

An ontology for user-system interaction modeling

Based on description logics, the ontology described above models the semantics underlying interactions between one user and his environment. As it is not possible to distinguish who activates a switch, we assume that the user stands alone. Obviously, our work should fit the need of a multi-user application, provided each log can be associated with an individual user id (e.g. services are requested from individual interfaces).

The ontology is made of six fundamental concepts: -User: someone who interacts with his environment using spread devices -Device: resource activated by user i) remotely or ii) from the same room -Basic Service: directly provided by both local and remote devices -System: interface to collect interactions composed of devices and services -Location: part of the environment that can contain devices and users -Effect: effect of services activation on the environment (e.g. luminosity)

This ontology is a useful support to achieve supervision 1 using as few sensors as possible. From the way the user interacts with his environment, we seek to infer his location and his intention related to the effect produced by services activation. Knowing where the user is and what is his intention, it is then possible to lead to service delivery. This is made possible by the formal semantic of description logics underlying ontologies, but also by rule-based reasoning capabilities.

1 By supervision, we mean not only monitoring, but also service delivery.

From a non-intrusive monitoring to service delivery

User tracking is an important context-awareness process, from which location can be inferred. In addition to ontology concepts described above, we add semantic relationships to link them. Figure 1 (left side) shows part of the ontology terminology (TBox ) for user location. The Door concept is defined with exactly two links isBetween towards two concepts Room. Given this logic representation, rules are expressed to infer properties such as isLocatedIn, isEntering and isLeaving following basic services activation. Consider for instance the rule 1 from table 1: a user ?u who is using a local activation device ?d which stand in room ?r, is himself in this room. Given the user location, it is then possible to deliver services based on the effect concept: when a user asks for a service, he knows that this activation will have an effect on his environment. Basic services are linked with effect (Luminosity, Darkness, Heat, Quietness, etc.) through the hasEffect property (cf figure 1, right side). The user intention can thus be inferred from this kind of property. Let's consider the rule 4 from table 1: the user ?u who is asking for a basic service ?b, wants the effect ?e associated with it.

Coupling the user intention with location information leads to service delivery. A user who wants luminosity (e.g. wants(User1, Luminosity)), and who enters a room will be offered services allowing to increase luminosity in the incoming room. Service delivery is based on suggestDeviceActivation and SuggestServiceActivation, applied respectively on System and Device (cf rules 5 and 6 from table 1). The design of such an ontology has to be flexible to adapt to user needs and evolving environments. On the prospect of a large scale service deployment, we choose not to do this complex modeling task by hand, but to automate it by means of a model-driven process.

A model-driven architecture for ontology design

In the context of ambient assisted living, we strongly believe that, in addition to the home automation expert, a non-expert -but someone who is able to properly consider the context of living with a disability (e.g. an occupational therapist, a family member) -must be included in the design process. We previously developped in [START_REF] Allègre | Model-Driven Flow for Assistive Home Automation System Design[END_REF] a model-driven flow for the design of assistive home automation systems. Platform independent modeling (PIM ) are separated from platform specific modeling (PSM )2 , so that non-expert designers can focus on specification rather than on implementation. Using a domain specific language, they can model both environment and interactions to define respectively the context of use and the user needs, before automatically generating the control code 3 .

In this paper, we propose to automatically generate an ontology with the same domain specific language. Figure 2 shows the supervision module that consists of several steps integrated in the existing design flow. In step 1', relevant model concepts defined by this latter at PIM level are transformed into ontological concepts (PSM level). They have some basic information such as the layout of housing, devices parameters, but also services defined to consider user needs. Step 2' is about generating a knowledge base conformed to OWL, the most commonly used ontology language. The joint use of a model-driven architecture and ontology4 makes sense to benefit from their respective advantages, thus leading to the definition of a complete control and supervision system for smart home. MDA and languages derived from this paradigm offer few means to formally interpret the semantic of models at runtime. On the contrary, ontologies The implementation of our ontology as a non-intrusive supervision model is illustrated as follows. When the user asks for a service, it is first processed by the home automation server, before sending him back its execution result. The service request is also translated into an assertion in the ontology that represents and shares context information coming from various home automation platforms. Then, an inference engine has to interpret this assertion and to enrich the knowledge base with i) monitoring information that an home automation non-expert will be able to analyze to define new rules or new services and ii) service delivery on the user interface based on information previously inferred.

Case study

To validate our approach, a person tracking application for monitoring and for dynamic service delivery is presented in this section. This case study is based on a real monitoring dataset recorded in the Domus smart home [START_REF] Kadouche | Support Vector Machines for Inhabitant Identification in Smart Houses[END_REF][START_REF] Chikhaoui | A New Algorithm Based on Sequential Pattern Mining for Person Identification in Ubiquitous Environments[END_REF] from the University of Sherbrooke (cf figure 3). It consists of data coming from sensors. Since our aim is to avoid the use of a maximum of sensors, we have converted some of these sensors logs by user requests on devices (cf table 2) when such devices can be remotely controlled 5 . We consider an excerpt of activities performed in day 2 by user 1 [START_REF]Domus Monitoring Dataset[END_REF].

Step 1: The user (user1 ) switches on the bedroom light: "[Log 1]: Li-ght5105.SwitchON " (log 1 from table 2). This log is introduced in the ontol- ogy through the following assertions6 : "[Assertion 1]: User1 isAskingForService SwitchON "and "[Assertion 2]: User1 isUsing Light5105 ". Knowing that the ontology defines Light5105 as a local activation resource located in Bedroom, and that this latter has a SwitchON basic service with Luminosity effect, rules 1 and 4 from table 1 can be applied. The inference engine execution allows then to lead to the following inferences: [Inference 1]: User1 isLocatedIn Bedroom and [Inference 2]: User1 wants Luminosity.

Step 2: The user now wants to go to the bathroom: after opening the bedroom door, he switches on the kitchen light: "[Log 3]: Light5101.SwitchON ". Since this light is considered as a local activation resource, the following inferences are realized: [Inference 3]: User1 isLocatedIn Kitchen and [Inference 4]: User1 wants Luminosity.

Step 3: The system proposes to switch off the bedroom light. This proposition is based on the inverse effect property (e.g. Luminosity isInverseOf Darkness). The user being located in the kitchen, the system proposes to switch off lights in other rooms. In this scenario, this service is accepted by the user: "[Log 4]: Light5105.SwitchOFF ".

Step 4: To go into the bathroom, he asks for the door opening: "[Log 6]: Door0216.Open". Considering this door linked with exactly two rooms with the isBetween property, rules 2 and 3 of table 1 can be applied. The user-door interaction allows the inference engine to extract the following information: [Inference 5]: User1 isEntering Bathroom and [Inference 6]: User1 isLeaving Kitchen.

Step 5: Finally, the system proposes him to switch on the bathroom light (log 7), from user intention and location previously inferred. These inferences are made possible trough rules 5 and 6 from table 1. Later, with the location information from sensor IR01, the system proposes him to switch off the bathroom light (as in step 3). This proposition is accepted by the user: "[Log 12]: Light5106.SwitchOFF ".

The Domus smart home is equipped with infrared movement detectors, pressure detectors, switches contacts, and flow meters. This case study shows that the automatic inference of the user location by our model should limit the use to 37.5% of the sensors that has to be deployed to ensure this person tracking application. Sensors that can be removed are essentially switches contacts replaced by user requests and infrared movement/pressure detectors replaced by logic inferences. This case study can be generalized to all situations where the user moves from room to room, and needs the same services throughout its path (e.g. listening to music, watching tv, etc.). In a context of severe disabilities without economic constraint, quality of services may be preferred. If we decide to keep all sensors from this case study, this approach can improve the system dependability by sensors information redundancy (66.7% in this case study 7 ).

Conclusion

The considerable increase of dependent people brings out important needs of monitoring to ensure a safe and secure independent living. In this context, this paper has presented a non-intrusive supervision model based on ontology for smart home. The supervision of such systems is not based on sensors use, but on the semantic analysis of logs recorded each time the user asks for a service. An ontology with associated rules allows the inference engine to extract information usually sensored. Finally, a case study based on a real monitoring data set validated our approach through a person tracking application.

This case study voluntarily uses very few sensors. A logic inference engine has to deduce the user location, thus replacing the use of movement detectors. However, to achieve a complete monitoring system over time, this approach should consider both user feedback and sensors whose information can not be inferred from home automation logs: information coming from biomedical sensors (e.g. blood pressure, core body temperature, etc.) or accurate information meters (e.g. temperature, luminosity, etc.). These latters could also be translated in assertions in the supervision model. Nevertheless, this approach allows the use of sensors to be limited, leading to a less intrusive and less expensive solution. Furthermore, it can improve the system dependability by redundancy of sensors information, especially important in ambient assisted living. In both cases, this approach provides a solution that may be more acceptable.

Service delivery is also an interesting contribution to allow users faster and easier access to services from their control interfaces 8 . A full-scale deployment is planned in two apartments of Kerpape MFRRC9 to evaluate how it can assist patients in their everyday life.

Fig. 1 .

 1 Fig. 1. Ontology (excerpts): user location (left side) and service delivery (right side).

Fig. 2 .

 2 Fig. 2. User-centered design flow for the control and the supervision of smart home

Fig. 3 .

 3 Fig.3. A person tracking scenario in the Domus smart home, from[START_REF]Domus Monitoring Dataset[END_REF].

Table 2 .

 2 Dataset (excerpt) of series 1 / user 1 / day 2 considering devices activation (*: logs coming from user requests /

	N •		Time	ID	Name	Location	Service/State
	1	*	8:56:17	5105	Light	Bedroom	SwitchON
	2	*	8:56:26	0215	Door	Bedroom	Open
	3	*	8:56:31	5101	Light	Kitchen	SwitchON
	4	*	8:56:35	5105	Light	Bedroom	SwitchOFF
	5	*	8:56:39	0215	Door	Bedroom	Close
	6	*	8:56:44	0216	Door	Bathroom	Open
	7	*	8:56:46	5106	Light	Bathroom	SwitchON
	8	•	8:56:51	Fl01	T.ColdWater	Bathroom	Open state
	9	•	8:56:54	Fl01	T.ColdWater	Bathroom	Close state
	10	•	8:57:09	IR01	Infrared	KitchenSink	Close state
	11	•	8:57:12	IR02	Infrared	KitchenOven	Close state
	12	*	8:57:12	5106	Light	Bathroom	SwitchOFF

• : logs coming from sensors).

In accordance to the Model-Driven Architecture (MDA) specification[START_REF] Miller | Model Driven Architecture (MDA). Object Management Group[END_REF].

Experiments have been conducted to validate the assumption made above.

The Ontology Definition Metamodel[START_REF] Consulting | Ontology Definition Metamodel[END_REF] is a specification to unify MDA and Ontology, making OWL conformed with the three modeling levels of the MDA specification.

As an example, Fl01 is a flow meter mounted on the cold water tap of the kitchen sink. It sends an event when the tap is opened but cannot be remotly controled

Each home automation log is divided into isAskingForService and isUsing properties.

On table 2, the information from

sensors over 12 can be replaced/supplemented by a service request (see logs with *).

Disabled persons with severe motion impairment (e.g. tetraplegia) can only control their home automation system by means of an automatic scrolling interface. In such restrictive situations, the selection of the desired service is very time consuming.

Kerpape Mutualistic Functional Reeducation and Rehabilitation Center.

⋆ This work is financed by the Brittany Region and is the result of a cooperation supported by Kerpape MFRRC (http://www.kerpape.mutualite56.fr/).