Automatic Extraction and Size Distribution of Landslides in Kurdistan Region , NE Iraq

This study aims to assess the localization and size distribution of landslides using automatic remote sensing techniques in (semi-) arid, non-vegetated, mountainous environments. The study area is located in the Kurdistan region (NE Iraq), within the Zagros orogenic belt, which is characterized by the High Folded Zone (HFZ), the Imbricated Zone and the Zagros Suture Zone (ZSZ). The available reference inventory includes 3,190 landslides mapped from sixty QuickBird scenes using manual delineation. The landslide types involve rock falls, translational slides and slumps, which occurred in different lithological units. Two hundred and ninety of these landslides lie within the ZSZ, representing a cumulated surface of 32 km. The HFZ implicates 2,900 landslides with an overall coverage of about 26 km. We first analyzed cumulative landslide number-size distributions using the inventory map. We then proposed a very simple and robust algorithm for automatic landslide extraction using specific band ratios selected upon the spectral signatures of bare surfaces as well as posteriori slope and the normalized difference vegetation index (NDVI) thresholds. The index is based on the contrast between landslides and their background, whereas the landslides have high reflections in the green and red bands. We applied the slope threshold map to remove low slope areas, which have high reflectance in red and green bands. The algorithm was able to detect ~96% of the recent landslides known from the reference inventory on a test site. The cumulative landslide number-size distribution of automatically extracted landslide is very similar to OPEN ACCESS Remote Sens. 2013, 5 2390 the one based on visual mapping. The automatic extraction is therefore adapted for the quantitative analysis of landslides and thus can contribute to the assessment of hazards in similar regions.


Introduction
Landslides are geological phenomena, which include a wide range of ground movements, such as rock falls, deep slope failures and shallow debris flows.These mass movements can happen in offshore, coastal and onshore environments.Although the action of gravity is the primary driving force for a landslide to occur, there are other contributing factors affecting the original slope stability [1].
Guzzetti et al. [2] considered that any mass movement is part of a landslide inventory map.In general, these maps show the spatial distribution of deposition and erosion areas produced by gravity-induced mass wasting which may vary in type, age and activity.In addition, landslide inventory maps involve landslide detection, recognition and classification.Size is one of the most important aspects in the recognition of landslides [3].The estimation of the size and probability of landslide occurrence is thus part of any landslide hazard assessment [4].Until recently, visual interpretation over aerial photographs and orbital high spatial resolution images remained the major source for landslide inventory map preparation.So far, visual interpretation of landslides with field checking is more accurate than automatic extraction in many cases.Nonetheless, the accuracy of automatic procedures has increased steadily and thus they are more and more adequate to be used to generate landslide inventory maps.
Borghuis et al. [5] detected landslides in southern Taiwan by using the Normalized Difference Vegetation Index (NDVI), supervised and unsupervised classifications as well as manual delineation.The authors indicated that unsupervised classification could detect 63.1% of the landslides, which were mapped manually.Nichol et al. [6] used multi-temporal images and showed that change detection is a successful technique to detect up to 70% of the landslides from SPOT multispectral data with spatial resolution of 20 m.The higher resolution images like pan-sharpened IKONOS imagery could be utilized for more detailed landslide inventory maps and the most affected areas.Martha et al. [7] used object-oriented methods (OOA) for semi-automatic landslide detection.In their investigation, they used high spatial resolution data and a Digital Elevation Model (DEM) to generate several parameters such as NDVI, slope, flow direction, hillshade, terrain, curvature, and stream network.They achieved a cognition accuracy of 76.4%.They then used the cumulative area-frequency distribution to compare manual delineation with automatic extraction of landslides.
A cumulative log landslide number distribution is a graphic plot that shows the relation between log size and log number of landslides.The x-axis represents log size and the y-axis represents log of cumulative number.The cumulative distribution allows quantification of the results of the inventory map by determining the statistics of landslide sizes.The cumulative distribution of landslides is given by (Equation (1)): (1) where: = cumulative number of landslides, α = the cumulative exponent, and = constant.Fujii [8] was probably the first researcher who investigated these types of structures.He studied 800 landslides in Japan caused by heavy rainfall [9][10][11].The cumulative number is calculated by using the area or a combination of area and volume of landslide triggers [12][13][14][15][16][17][18][19][20].Other studies have been made using only the volume of landslides [9,21,22].
In this paper, we explore the potential of very high spatial resolution data and a semi-deterministic spectral characterization of scars in order to increase the efficiency of automatic extraction.First, we present the study area, the available remotely sensed data and the data processing.Second, we discuss the advantages and limitations of the proposed method for the recognition and mapping of landslides.Finally, we present and then discuss the results of frequency-size analysis for two zones of the Zagros orogenic belt.
This study involves four main steps of data processing.Firstly, we prepared a landslide inventory map of delineated landslides manually from QuickBird images.Secondly, we generated cumulative frequency-size distribution formulas of two tectonic zones of the Zagros orogenic belt.Then, we improved the automatic detection of landslides in the region.Finally, we made a statistical comparison between the automated extraction and the manual delineation for parts of the study area.

Study Area
The study area lies within the Kurdistan region in the northeastern part of Iraq (36°30′N-35°30′N, 44°26′E-46°21′E).It comprises the Iraq Zagros Mountains, in which mass movements threaten many towns and big villages.The study area covers about 6,554 km² (Figure 1).

Geological Setting
The NW-SE trending Zagros orogenic belt is approximately 2,000 km long, extending from southeastern Turkey through Iraq to southern Iran and part of the Alpine-Himalayan mountain chain [24,25].The Iraqi territory is divided into four main tectonic zones of different characteristics (rock types, age, thickness and structural evolution).These four zones are (1) the Inner Platform (stable shelf); (2) the Outer Platform (unstable shelf), which includes the Imbricated Zone (IZ), the High Folded Zone (HFZ), the Foothill Zone and the Mesopotamia Foredeep; (3) the Shalair Zone, and (4) the Zagros Suture Zone (ZSZ) (Figure 1) [23,[26][27][28][29][30][31].The Suture Zone is characterized by ophiolites [26].The study area consists of ten units formed during the Late Cretaceous and Mio-Pliocene periods due to the obduction and collision of the Arabian with the Iranian plate.Three NW tectonic zones formed in the ZSZ: The Qulqula-Khwarkurk Zone, which consists of radiolarian cherts, mudstone, and limestone, conglomerates and basic volcanic rocks of increasing occurrence to the NE of the study area; the Penjween-Walash Zone, which consists of metamorphosed rocks, carbonate beds with volcanic and pyroclastic rocks; and the Shalair Zone, which consists of meta-pelitic and meta-carbonates, volcanic, and metamorphosed rocks.These structures are overlain by different types of Quaternary sediments (river terraces, alluvial fans, slope sediments, valley fills, flood plains, polygenetic sediments) [32,33].The HFZ and the IZ include sedimentary rocks, mainly limestone, but also dolomite, sandstone, claystone, shale, conglomerate, marl, siltstone and bitumen, which range from Triassic to Quaternary [34].
The main thrust fault forms the NE to NW tectonic contact between the IZ and HFZ [26].Due to the similarities of the IZ (a ~25-km-wide narrow belt) with the HFZ [34], we considered it as a part of the HFZ.On the northwestern side, the HFZ includes more than 42 mainly NW-SE oriented, asymmetrical, anticlines and synclines.Most folds are box shaped with steep SW or W flanks.They were influenced by transversal fault and reverse-slip fault systems.On the other side, within the Shaliar Zone in the southeastern part of the study area, the ZSZ includes few E-W or NE-SW oriented anticlines and synclines.
The study area is characterized by rugged topography in most parts.The elevations range from 376 m to 3,415 m with slopes ranging between flat and 81°.This landscape is the result of the combined effects of tectonic uplift, erosion, and different rock strengths.The actual geomorphological units developed under severe erosion and moderate weathering [32,33].Most of the rocks outcropping in the region are characterized by a high brightness after they have been freshly exposed.

Climate
The Lesser Zab River collects all the surface water that enters the basin as precipitation (rain and snowfall) and flows to the Dukan Lake.The Köppen-Geiger climate classification system [35] characterizes the climate as warm temperate with dry and hot summer (Csa).The climate of this mountainous basin is defined by large seasonal variations in precipitation, temperature and evaporation, represented by dry summers and wet winters.The climatic data is gathered from the Agro-Meteorological Department of the General Directorate of Research & Agricultural Extension of the Ministry of Agriculture of the Kurdistan Regional Government).Most of the entire annual precipitation (895.3 mm) occurs during October to May.The highest average monthly precipitation value of the last 10 years was measured in January with 198.3 mm.Monthly average temperatures vary between 0 °C (February) to 38.4 °C (August).Above 1,500 m, heavy snowfall occurs in winter.Heavy rainfall, heavy snowfall and rapid snow melting subsequent to sudden change in temperature, lead to the incidence of landslides in spring.In this area, the annual evaporation is up to 1,939.7 mm.The highest average monthly evaporation of the last 10 years was measured in July with 368.1 mm (Figure 2).

Landslides
The study area is susceptible to frequent landslide occurrences due to both natural and human-induced reasons.The main natural factor, which causes the landslides in the study area, is the very rugged topography across the elongated anticlines.Other natural factors are: the rapid snow melting in spring, the focused precipitation, and the relatively heterogeneous geology and geomorphology.Human-induced causes of slope failures include the construction of the Dukan Lake dam, civil engineering activities like road cuts, overloading of the top or undercutting of the toe of slopes and unsuitable agricultural practices.
Landslides and other types of mass movement are related to water supersaturation, which increases the internal water pore pressure.Among the numerous events that occurred in the region, a massive and good example is the landslide that occurred 10 km eastwards of Chowarta town, within the Gora Kachaw Mountain south of the study area [36].It damaged parts of the village Galley and led to the death of 20-25 persons.It was caused by high water saturation of the ground due to ice melting.The heavy snow and rain caused mud banks to collapse, and the mudflow buried houses and blocked roads.On 24 February 2013 part of a cliff in Sara Mountain collapsed due to the heavy rain and blocked the road from Dukan to Bangard towns.
Along the main road to Gimo Mountain, after Kanaro village, large areas are affected by landslides.A few hundred thousand cubic meters of igneous rock of the Mawat massif has collapsed in the last decade.Rock falls and toppling has occurred in different lithological units along steep slopes and gorges independently or collectively.Recent events caused road blockings and several towns nearby are threatened and regularly affected-a large rock fall was witnessed recently in the Chowarta town along the road to Gimo village near Dukan Lake (Figures 3 and 4   Translational slides or rockslides are common in hard rock and occur along planar shear surfaces [37].This type of slide damages more than rock fall.Slumps or spoon-like landslides are common in weak layers, especially when they have been softened by percolation of rainwater [37].

Data, Data Preparation and Workflow
The availability of new satellite sensors with better spatial and spectral resolutions made the use of satellite data more attractive than aerial-photographs to detect and investigate landslides.For this study, 60 cloud-free QuickBird scenes were used.They were obtained via the Ministry of Planning (Iraq) and were acquired from 18 August 2002 to 29 August 2006.The final product is a 3-bands, 8-bit, 0.6 m resolution mosaic covering more than 6,554 km 2 .It has three visible spectral bands in the blue (450 to 520 nm), green (520 to 600 nm) and red (630 to 690 nm) and one panchromatic band (450 to 900 nm).The mosaic did not include the near infrared band (760 to 900 nm).The three visible spectral bands data were pan-sharpened using the UNB algorithm [38].The mosaic was radiometrically corrected, ortho-rectified and projected using theWGS84 data and the UTM 38N projection.In addition, we used eight ortho-rectified and radiometrically corrected Advanced Spaceborne Thermal Emission and Reflection (ASTER) and their extracted stereogrammetric digital elevation models (DEMs) ASTER using Nadir (N) and backward looking (3B) red bands.The ASTER level 1A scenes were acquired between 24 August 2003 and 15 September 2006 with a resolution of 15 m.The ASTER data was necessary because of the unavailability of the NIR band at higher spatial resolution.
Twenty previous geological maps from published reports of [39][40][41][42] were used in this study (Figure 6).These maps have varying scales and cover different parts of the study area.Bolton [39] provided a geological map of 1:100,000-scale that covers the northeast part of the study area.In this report, 15 major landslides were detected from 1:20,000-and 1:40,000-scale aerial photographs and additional fieldwork.Paver et al. [40] produced 15 geological maps at a scale of 1:20,000 covering the east and northeast of the study area.The report lists 25 major landslides based on fieldwork observations.Buday [41] delivered four geological and geomorphological maps at a scale of 1:100,000, which cover the east part of the study area.The report lists 39 major landslides based on fieldwork observations.Abdulaziz [42] published six engineering geological and geomorphological maps at a scale of 1:100,000, enclosing the entire study area.They delineated 117 major landslides based on 1:40,000-and 1:50,000-scale aerial photographs, and fieldwork for part of the study area, and a hardcopy of MSS Landsat imagery (1:250,000) for the rest.Environment for Visualizing Images (ENVI) software was used to perform the data operations (mosaic, subset, threshold, and automatic landslide extraction).All GIS operations (base map preparation, raster vector conversion and area calculation) were done using ArcGIS10 and all statistical operations were done using R software.

Methodology
The manual delineation of landslides using two sources provided the inventory map (Figure 7).Firstly, previous geological maps [39][40][41][42] were scanned with 300 dpi and georeferenced to the UTM coordinate system in Zone 38 north.Then, the landslides were digitized in a GIS.Secondly, visual interpretations and digitizing of QuickBird scenes (partly verified by field survey in different parts of the study area) were carried out.The landslide boundaries were identified based on attributes such as tone, texture, headwall scarps, and associations like the pathway of the material movement and fragments of transferred materials from the satellite images.The two databases of digitized landslides were matched to create the final inventory map.

Automatic Landslide Detection
We developed a robust methodology to analyze and evaluate the landslides in the study area (Equation ( 2)).This methodology is neither meant to replace standards methods nor to be universal.On the other hand, the algorithm is easily adaptable to other (semi-) arid areas (like our study area; Section 2.2).The idea of the brightness indicator is based on the relatively high reflectance contrast between landslides and their background.Particularly, the high reflectances of landslides in green and red bands allow their detection [6].The rock types are composed of carbonate rocks (limestone and dolomite), marl, claystone, shale, sandstone, conglomerate and chert in this area.Almost all of these rocks reveal a higher reflectance in red and green bands after they have been recently exposed.

Brightness indicator 255
Green and red bands are stored as the second and third band in QuickBird data.The brightness indicator as described in Equation ( 2) can be described as a brightness of both red and green bands.Additional portable Analytical Spectral Devices (ASD; [43]) FieldSpec measurements were also performed in order to support the discussion of the results obtained from the remotely sensed images.The samples were collected from different landslides in the study area.The samples encompass seven selected rock samples (i.e., four limestones, two conglomerates and one of chert).The measurements were conducted in the Iraq geological survey laboratory.
Figure 8 shows the high reflection of carbonate minerals (calcite and dolomite), clay mineral (muscovite, montmorillonite, and kaolinite) and quartz with a wavelength of 0.45-0.7 µm.The outcrops of rocks due to erosion reduced the reflectance in red and green bands because of the oxidation.In addition, the moisture content further reduced the reflectance because of the absorption of water in green and red bands [44].We used a test area on which we have a complete inventory.The satellite data of the selected area was acquired free of any snow, cloud, salt and sand dune coverage.The brightness indicator was developed to automatically extract and generate the inventory map for recent landslide areas.The brightness indicator is a raster map with values ranging from zero to 255 (Equation ( 2 Houses, unpaved roads, and salty areas are usually characterized by low and flat slope.Therefore, we apply a slope threshold.The selection of the threshold values for both brightness indicator and slope data was determined experimentally by using a contingency table [46].Additionally a further mask was applied in order to remove the vegetated area (i.e., yellow plants, and grasses) using the normalized difference vegetation index (NDVI).The NDVI is expressed as the difference between the near infrared (NIR) and red bands normalized by the sum of those bands (Equation (3), [47]).NDVI is a common tool to detect the landslides in distinct environments [5,7]. (3) We calculated NDVI using both red and NIR spectral bands from ASTER.This procedure was necessary because the NIR band was precluded from the QuickBird dataset.The final thresholded brightness indicator raster was then used to generate the areas representing potential landslide surfaces.

Statistics
After calculating the landslide size, histogram distributions were plotted to show the frequency of occurrences throughout the study area.Categorization was made based on tectonic zones in the HFZ and the ZSZ.The cumulative distribution number of landslides was calculated using 30 bins.The relation between log size and log number as well as the density of the log of landslide size were plotted.The kernel density was adopted to estimate the probability density of log size.The cumulative distribution number of landslides >100 m 2 often follows a power law between the cumulative number of landslides (N Cl ) and the landslide area (A) [8].The visual delineation was used as a reference to validate the automatic extraction model.We used a contingency table and the Equations ( 4) and ( 5) ( [46]) to determine the accuracy of the automatic extraction.In addition, the contingency table was employed to select the best values of the slope angle and the brightness indicator values, located in non-vegetated area.Both types of landslides extracted (manually and automatically) are delineated to one element of the set of positives (there is a landslide) and negatives (there is no landslide) (Figure 9).[46]).

Automatic
False Negative True Negative From the input vector layers (manual and automatic extracted landslides polygons), there are four possible results.If the landslide manually extracted is positive and the landslide automatically extracted is positive, it is considered a true positive; if the landslide manually extracted is positive and it the landslide automatically extracted is negative, it is considered a false negative [46].If the landslide manually extracted is negative and the landslide automatically extracted is negative, it is considered a true negative; if the landslide manually extracted is negative and the landslide automatically extracted is positive, it is considered a false positive [46] (Figure 12(e,f)).The contingency table was performed using 'Boolean operations' between the polygons detected manually and automatically.

Positives correctly classified
Total positives (4) (5) The positive correctly classified is true positive; the negative incorrectly classified is false negative.The total positives are the sum of true positives and false negatives.The total negatives are the sum of the false positives and the true negatives [46].An error table is a good way to detect the error in both types: commission (error of inclusion) and omission (error of exclusion) [48].
In addition to a contingency table, two tests were used to determine the accuracy of the automatic extraction.The first is expressed in percentage and represents accuracy between the number of overlapped landslides and of manual delineation (Equation ( 6)).The calculation was made from Equation ( 6), where the overlapped landslide number represents the intersection of landslides from the manual and the automated method.accuracy assessment overlapped landslide number manual delineation landslide number 100 Secondly, we plotted the relation between log size and log number for the automatic and the manual landslide delineation.An exponent of both was calculated using R software to show the accuracy of automatic landslide extraction.

Manual Delineation
Since the estimate of landslide size is part of landslide hazard assessment [4], we digitized 196 large landslides from twenty previous geological maps [39][40][41][42] (Figure 6).From the QuickBird images 3,190 landslides were derived using visual interpretation (Figure 3); 290 of them were in the ZSZ and 2,900 landslides were in the HFZ.The total landslides coverage accounted for an area of 58.5 km 2 .The detected landslide size varied from 100 m 2 to 3.33 km 2 .The density of landslides per km 2 for the entire area was 0.47 (Table 1).The majority of the landslides occur in the northern part (Figure 3).The probability density distribution of 290 landslides within the ZSZ and 2,900 landslides within the HFZ are plotted in Figure 10(b).The bandwidth estimation of kernel density for the HFZ is 0.14, and reflects large variance in landslide sizes.For the ZSZ the kernel density is 0.28 and indicates that landslide size variance is smaller compared to the HFZ. Figure 10b shows that there are a limited number of very large landslides in the HFZ.On the other hand, there are a limited number of very small and very large landslides in the ZSZ.The maximum distribution of landslides within the ZSZ ranged from 0.001 to 0.1 km 2 , whereas the landslides within the HFZ reached from 100 to 10,000 m 2 .
The largest 10 landslides (>1 km 2 ) are located in ZSZ, within Naopurdan-Walash Group (Paleocene-Oligocene).Naopurdan-Walash Group is a composite of gray shale with thin beds of green greywacke, and lenticular volcanic rocks.These 10 landslides (largest red patches in Figure 3) are the largest landslides in the study area covering a total area of 20.13 km 2 , located in the northeast-east of the study area.The fit of a power law to the cumulative distribution of event sizes is commonly used to analyze landslide [12][13][14][15][16][17][18][19][20].Here, we calculated the cumulative number-area of the landslides within the HFZ and the ZSZ.The cumulative distributions of 290 landslides within ZSZ and 2,900 landslides within HFZ show statistical correlations >0.94 (Figure 11).The power exponents were α = −0.3 for the ZSZ, and α = −0.9 for the HFZ, for landslide sizes <1,000 m 2 .However, the exponent of landslides (>1,000 m 2 ) were α = −0.68 for the ZSZ, and α = −1.1 for the HFZ.

Automatic Landslides Detection and Validation
We created four thresholds of brightness indicator values (84, 86, 110 and 160) and four thresholds of slopes (5, 10, 12 and 15) (Table 2).We created 10 vector layers devoid of vegetation and displaying a representative variation of the thresholds (Table 2).These 10 maps were filtered from vegetation using the ASTER derived NDVI map.The map which has non-vegetated areas, slope >10° and the brightness indicator >84 gives the best results for the automatic extraction (Table 2).For this map, the area of the true positives (intersect between automatic landslide and manual landslides detected) was ~1.4 km 2 , representing ~83% of the visually interpreted ones.The validation of the brightness indicator >84 and slope >10° non-vegetated map indicated more errors of commission rather than errors of omission.The omission in the true landslide areas for this map was ~17%, and the commission was ~31%.The examination of the visual interpretation proved that areas, erroneously commissioned by the automated methods, were found to be bare soil, weathered bedrock or quarries.Of recent landslides 1,798 were derived using manual extraction.The brightness indicator was able to detect 1,726 of these landslides, which represent an accuracy of 96%. Figure 13 shows the distribution plot of automatic extraction and visual interpretation of recent landslides within the selected area.The two lines show a statistical correlation superior to 0.92.The exponent in the plot of Figure 13 reveals α = −0.63 and α = −0.58 for the automatic extraction and the manual delineation, respectively.The distribution matching between the automatic and manual exponent was 92.7%.

Discussion
The cumulative distribution plot (Figure 10) shows that the landslides in the HFZ are less variable in size than in the ZSZ.This is due to the difference between rock types and tectonic activities.The landslide inventory map shows that larger landslides are older and difficult to map automatically.Old landslides allowed the reestablishment of vegetation.
Visual interpretation of remotely sensed data for landslide detection is more time consuming than automatic procedures.In addition, visual interpretation is based on the interpreter's perception and experience to detect landslides.It is difficult to identify small landslides with the manual delineation of QuickBird data at a spatial resolution of 0.6 m.The brightness indicator was able to detect smaller landslides than the manual delineation.
The ASD derived reflectance in the 0.45-0.7 µm range showed significant differences between the fresh and weathered faces of the three rock types (i.e., limestone, conglomerate and chert).In a similar way, fresh faces of the three rock types showed similar reflectance to the recent landslides detected by the QuickBird data.Weathered faces showed similar reflectance to the bed rock, which surrounds the recent landslides.The conglomerate and the limestone fresh face showed higher reflectance than the weathered face.The contrast between recent landslides and surrounding area corroborates with previous investigation conducted by Nichol et al. [6].However, this assumption is only valid for most of the recent landslides, which consist of the rocks that have higher reflectance for fresh face than the weathered face in green and red bands.On the other hand, the chert fresh face showed lower reflectance than the weathered one.Therefore, we cannot always consider that the high reflectance of both green and red bands of the recent landslides allow their detection as reported by [6].
The threshold values of brightness indicator values and slope can be different from place to place but their selection is easy to retrieve experimentally.The algorithm is robust and able to recognize the recent landslides.The brightness indicator values used here could be applied to monitor landslides in other regions, with similar characteristics (geology, land cover, etc.).We cannot consider this technique to be a general solution to detect landslides but the modified (semi-automatic) algorithm could help to produce a primary inventory map especially after earthquakes.Some errors in the method are due to the (i) the difference in spatial resolution between ASTER and QuickBird; (ii) difficulty to extract vegetation from the landslide polygons (false positive area); (iii) ambiguity between roads and tilted roof with slopes >10°) landslides; (iv) landslide occurrences in dark rock (like chert; true negative).The trace of landslides in shadow areas and wet lands is difficult to distinguish because of low brightness indicator values (true negative).Land surface with false high values (e.g., high spectral response from houses, salt, snow, cloud coverage and both unpaved and cemented roads; false positive) can also be confused with landslides.Masks are a good solution to overcome the confusion of the land surface with false high values.Due to re-vegetation of old landslides, the brightness indicator values do not allow detection.
From the total detection of landslides with the automatic extraction ~96% was made by visual interpretation.The accuracy assessment between the manual delineation and the automatic extraction was 92.7% (Figure 12).The difference between cumulative landslide number-size distributions of the automatically and visual landslides extracted was either due to the separation of single landslides, by one segment of the manual delineation shapefile or by two or more segments of the automatic extraction shapefile (Figure 12(e)) (ii) merging of two or more landslides, represented by two or more segments of the manual delineation shape file, to only one segment of the automatic extraction shape file (Figure 12(e)).
Landslide mapping by field investigations has been shown to be a challenging task in mountainous terrain.The remote sensing techniques and the high spatial resolution of the QuickBird and ASTER images employed in this study provide a viable tool for landslide detection.

Conclusions
Over 3,190 landslides including different types of landslides (rock falls, debris flow, translational slides and slumps) have been mapped in this study by using QuickBird images.Among these, most of the landslides occurred in the northern part of the study area.The Naopurdan-Walash Group within the Zagros Suture Zone (ZSZ) involved the largest 10 landslides.The exponent of landslide size distribution is <1,000 m 2 with α = −0.3 and α = −0.9 for the ZSZ and the HFZ, respectively.However, landslides >1,000 m 2 have exponents (α = −0.68,α = −1.1)which means that the landslides in the High Folded Zone (HFZ) have a lower size than in the ZSZ.
The brightness indicator was applied in parts of the study area.This map was filtered by slope and NDVI maps to remove the low slope and vegetated area.The contingency parameter shows that the best map for the automatic extraction has non-vegetated areas, slope >10° and the brightness indicator >84.This map was detected ~1.4 km 2 .The proposed method detected about ~96% of the total number of recent landslides crosschecked by visual interpretation.The landslide number-size distribution analysis shows that the exponent was α = −0.63 and α = −0.58 for the automatic extraction and manual delineation, respectively.The error in the distribution matching was 7.3%.This was due to the different spatial resolution between used data and also to the difficulty to exclude vegetation, which represent only a few areas within some of the polygons of landslide.In addition, few types of land surface with high DN values show false positive landslides.Furthermore, the manual delineation was based on the level of expertise of the interpreter and the accuracy of digitizing making it vulnerable to human error.Nonetheless, the robustness and the facility to adapt the algorithm to varied regions make it a valuable tool for rapid hazard assessment.

Figure 2 .
Figure 2. Monthly precipitation, evaporation, maximum and minimum temperature, in the study area based on data from 2000 to 2006 (the Agro-Meteorological Department of the General Directorate of Research and Agricultural Extension of the Ministry of Agriculture of the Kurdistan Regional Government).
(S1,S2)).Any reactivation in this cliff will certainly destroy the town.Big blocks of heavily fractured limestone fall from a steep cliff near Dukan Lake (Figure4(S2)).The volume of these blocks of limestone reaches >1,300 m 3 .

Figure 3 .
Figure 3. Landslides inventory map of 3,090 landslides identified through visual interpretation of 60 QuickBird scenes (the S1-S9 markings correspond to sites mentioned in Section 2.3 and Figures 4 and 5).

Figure 4 (
S3) shows the translational sliding which threatens the new Dukan town-Sulimanyah road cutting the northeastern flank of the Surdash Anticline near Dukan Lake (in Figure 3 (S1)).The roads from Sulimanyah to Chowarta and Mawat are also affected by recurrent landslides of different sizes.

Figure 4 .
Figure 4. Typical examples of landslides within the study area.(a,b) Rock fall in station S1 and S2; (c,d) Rock slide in the study area near Dukan Lake at station S3.

Figure 5 (
S4-S6,S9) shows that slump sliding threaten roads in different sites of the study area.Clastic debris flows are also common in the weak layers.Especially, in the north-northwest of the study area, the debris is usually fine-grained.Numerous landslides develop NE of QalaDiza town.One of the largest one is still active and covers an area of 1.4 km 2 (Figure5(S5)).

Figure 5 .
Figure 5. Red, Green, Blue (RGB) 321 of QuickBird data showing examples of detected rotational landslides in the study area.The four images represent the landslides S4, S5, S6and S9, which were checked in the fieldwork (see Figure3).

Figure 6 .
Figure 6.The inventory map of the 196 major landslides digitized from previous geological projects [39-42].
The spectral range of the digital number (DN) of brightness indicator is between 0 and 255.The final products QuickBird data has only 8-bit.Therefore, the constant number (255) in the brightness indicator equation represents the radiometric resolution of QuickBird data (8 bit).
)).Low brightness indicator values indicate areas without landslides (e.g., vegetation, rivers, and paved roads).On the other hand, high brightness indicator values indicate recent landslides, but also roofs of buildings made of clay and lime or just clay, cemented roads covered by clay and lime and unpaved roads, covered by clay, yellow plants and grasses.

Figure 10 .
Figure 10.Statistic plots, (a) Histograms showing the number (count) and size of landslides in logarithmic coordinates; (b) Probability density distribution showing 290 landslides within the Zagros Suture Zone (ZSZ), and 2,900 landslides within the High Folded Zone (HFZ).

Figure 12 (
Figure12(a,b) shows examples of extracted landslides.Figure12(c,d) shows the brightness indicator values of these landslides (gray to light gray color), whereas Figure12(e,f) shows the manual delineation and the brightness indicator >84, which has slope >10°.

Figure 12 .
Figure 12.Typical examples of automatic landslide extraction: (a,b) landslides in QuickBird image represent S8 and S7, respectively (c,d) brightness indicator (e,f) relation between shape file, which extracted from the brightness indicator with DN number >84 and slope >10° and manual extraction of landslide.

Figure 13 .
Figure 13.Distribution plot of automatic extraction and visual interpretation of recent landslides within the selected area.In the graphs, the x-axis shows the log of size, the y-axis shows the log of the number of landslide.

Table 1 .
Statistical characteristics of the landslides (rotational landslide, flow, rock fall and rockslide).

Table 2 .
Contingency table parameters and rates for different slopes and brightness indicator non-vegetated maps.