Carmille Marie J. Escapea,b,*
aNationwide Operation Assessment of Hazards, U.P. NIGS, C.P. Garcia Ave., U.P. Diliman, Quezon City, 1101 Philippines
bNational Institute of Geological Sciences, University of the Philippines, C.P. Garcia corner Velasquez street, U.P. Diliman, Quezon City. 1101 Philippines
*Corresponding author: Email address: [email protected]
As Typhoon Koppu traversed the northeastern provinces of Luzon last October 2015, massive debris flow events were triggered in the foothills of the Sierra Madre range, causing devastation in the municipalities of Gabaldon, Laur and Bongabon in Nueva Ecija. These fast moving hazards had carried boulders as large as 6.5 m in diameter and had buried houses and destroyed structures situated within the alluvial fans where human settlements had thrived. Fortunately, early warnings had been issued and the citizens were able to act accordingly, thus no casualties were reported. Though it was a disaster that had been averted, this event could happen again and thus, it is important to study the deposits and their extent in order to gain further data that could be useful in scenario based hazard modelling. The extent and character of the debris flow deposits were identified from processing pre and post-typhoon Landsat 8 images. Delineation of the deposit was done through change detection via image differencing and thresholding applied to NDVI indices and Tasseled Cap Transformation indices. The differentiation between boulder to pebble and finer sediment deposits, as well as stream delineation was done by unsupervised image classification of the NDVI index of the post-typhoon image and the 3rd principal component produced from the Principal Component Analysis of the stacked red and NIR bands of the pre and post-typhoon images. Accuracy assessment based from field data showed that the methodology was able to give a good delineation between the different types of deposits.
Last 14 October 2015, Typhoon Koppu (Local name: Lando) entered the Philippine Area of Responsibility as a tropical storm, bringing with it a maximum sustained winds of 65 kph . By 17 October, it had escalated into a Category 4 Typhoon (based on the Saffir-Simpson hurricane wind scale) with maximum sustained winds of 240 kph . On the afternoon of the same day, the super typhoon center had made landfall over Casiguran, Aurora where its maximum sustained winds had reached 185 kph with central pressure of 920 hPa . As it hit the Eastern Luzon provinces, it had brought heavy to extreme rainfall that necessitated the raising of Public Storm Warning Signals 2 – 4 within that region . By the time it had weakened into a Low Pressure Area as it moved through Balintang Channel on the afternoon of 21 October, the typhoon had already caused severe flooding, numerous landslides including a huge debris flow event in Nueva Ecija, and a slew of casualties. It had caused great devastation in infrastructure and agriculture in Regions I-V and CAR (Cordillera Administrative Region), equivalent to 10 billion pesos. Over a 130 thousand houses were damaged and 48 people had been reported dead, 83 injured and 4 missing . Because of its calamitous effect, the Philippine Atmospheric, Geophysical and Astronomical Services Administration (PAGASA) had decided to decommission the name .
1.1. Geomorphology of debris flows and alluvial fans in Nueva Ecija
At the foothills of the Sierra Madre mountain range, the locals within the municipalities of Gabaldon and Laur, Nueva Ecija had reported experiencing flooding that later on brought a mixture of sediments, boulders and other detritus. Fortunately, they were able to immediately evacuate and there were no recorded casualties . This rapidly moving mass of slurry, soil, rock, water, air and other components (e.g. uprooted trees) is a geological hazard known as a debris flow. They are typically the result of the mobilization of poorly consolidated materials situated at the elevated regions of a drainage area. The motility could either be caused by rapid snowmelt runoff or torrential rainfall .
This debris flow event which affected the municipalities of Gabaldon, Laur and Bongabon in Nueva Ecija, was deposited in a geomorphic feature known as an alluvial fan. These fan-shaped features are typically present in valleys contiguous to a mountainous region where a channel emerges. They form by the accumulation of sediments that are deposited due to the rapid loss of carrying capacity that a channel or transporting medium experiences once it stems out from the highlands . When two or more neighboring alluvial fans converge, they are collectively called as a bajada . Because they occur in piedmont zones, have channels that serve as water source and often have long periods of stasis or inactivity that could span from decades to centuries, alluvial fans had become ideal location for urban settlement. Unfortunately, they are also sites for two very disastrous hazards – alluvial fan flooding and debris flow, which often takes people by surprise because of the seemingly benign appearance of the topography .
Due to the presence of the Philippine Fault that transects the provinces of Aurora and Nueva Ecija in a northwest-southeast trend, the geomorphology of the affected region was conducive for the formation of alluvial fans and bajadas. From the mountainous region of Sierra Madre, numerous creeks emerged into the flatland traversed and formed by the fault. The rapid change in elevation resulting to progressive loss of flow power, availability of sediment source from high relief and short transport distance, allowed for the formation of alluvial fans and bajadas for almost each of the creeks (Fig.1(A)).
1.2. Debris flows in Nueva Ecija
As Typhoon Koppu raged into eastern Luzon, it had induced numerous landslides in the Sierra Madre range. The loose deposits of these landslides became the source material for the debris flow (Fig.1(B)). The torrential rainfall had also oversaturated the soil with water, allowing it to be easily mobilized at great velocities. Debris flows occurred at almost all of the alluvial fans situated in the mentioned municipalities. According to interviews conducted during a post-disaster fieldwork (28-30 October 2015), the debris flow events occurred between 9 AM to 12 Noon of 18 October. Boulders as large as 6.5 m in diameter were observed and a maximum thickness of 3.5 m of debris flow deposit was measured in Barangay Ligaya, Gabaldon (Fig.1).
Figure 1 illustrates the extent of the debris flow events. Although the deposits of debris flows are usually poorly sorted, the distribution of coarse-grained (boulders to pebbles) sediments versus that of finer grained sediments in the proximal and distal part of an alluvial fan can give information about the mechanism of flow  within the event . The flow mechanism can explain the various characteristics of a debris flow deposit such as its geometry, presence of inverse grading, imbrication and friable sediment blocks and clumps. It is therefore important to be able to characterize the debris flow deposits in order to understand the behavior of an event, upon which the observations gathered could be used as parameter for modelling scenario-based debris flow hazards.
It was therefore the aim of this study to identify and delineate various debris flow deposits through different remote sensing techniques applied in pre and post Typhoon Koppu Landsat 8 satellite images.
Various techniques such as object based image analysis, use of different indices, change vector analysis and image differencing, had been employed to interpret landslides from satellite images . Since debris flows are also considered as mass wasting events like landslides, these techniques can be used to interpret them. In fact, numerous studies regarding interpretation of landslides using optical remote sensing also included debris flows . However, characterization of the debris flow deposits using remote sensing techniques had not been given much focus. Therefore, the methods (Fig.3) used in this study were fused together from prior studies that aimed to characterize landslides  .
The satellite images used for this study were Landsat 8 images with spatial resolution of 30 m for the multispectral bands and 15 m for the panchromatic band. It has a 16 day revisit capability and can be freely downloaded from GloVis, EarthExplorer and LandsatLook Viewer. Three images were used; two of them served as pre-typhoon images (acquisition date: 06 September 2015 and 08 October 2015) while the other one served as post-typhoon image (acquisition date: 24 October 2015) (Fig.4). The acquired images are Level 1T images which are terrain corrected and are in the form of quantized, calibrated and scaled Digital Numbers (DN).
2.2. Data Preprocessing
In order to use the images, they were first transformed from Digital Number to Top of the Atmosphere reflectance using the equation provided by USGS (U.S. Geological Survey) . These are the reflectance values with contributions from clouds and atmospheric aerosols and gases. Atmospheric correction was done by the application of Dark Object Subtraction (DOS)  which built upon the assumption that there are pixels within the image that are under a complete shadow, and that their radiance received by the satellite is due to atmospheric scattering (path radiance) . After the correction, cloud masking was performed using the Quality Assessment (QA) Band which has pixels that contains integers that reflect surface, atmosphere and sensor conditions. Using the table provided by USGS and shown in Fig.1, cloud masking was performed automatically. For the pixels with uncertainties regarding presence of clouds, it was necessary to manually review whether they are indeed clouds or not.2.3. Change detection techniques
For change detection, image differencing was the primary method that was used. However, the Principal Component Analysis (PCA) which is a multi-temporal visual composite change detection method was also employed, but its product was used for image classification. PCA is a geometric transformation technique that is used for data reduction. A set of correlated variables are linearly transformed into a set of uncorrelated variable (principal components) in a different orthogonal coordinate system  . It compresses the spectral components by identifying components from which the most variance can be attributed. PCA was applied by stacking the images of pre and post red and NIR bands for a total of four components. Visual inspection showed that the 3rd component depicted the most variance within the debris flow deposits. Due to this variance, it was the principal component (PC) used for image classification.
For the image differentiation, another data transformation was used. The Tasseled Cap Transformation (TCT)  is similar to PCA in such a way that it takes linear combinations of the bands and summarizes and interprets them into a new set of bands. However unlike PCA, TCT’s first component is an index of the overall brightness of the image. The second component is for the greeness while the third component is an index of wetness. The rest of the bands usually contain “noise” and are not as useful. The TCT brightness, wetness and greeness index were produced for the pre(06 Sept) and post(24 Oct) satellite images. It was done by calculation using the equation derived by Baig (2014) and was applied to bands 2-7 (Eq.1) (Fig. 5).
Another index that was produced is the Normalized Differential Vegetation Index (NDVI) which measures the health and estimates green biomasses . NDVI is very helpful in distinguishing vegetation from bare earth and had been widely used for landslide detection and change analysis . It is calculated by using the red and NIR bands of the satellite image (Eq.2).
In order to show the change between the pre and post typhoon images, image differencing was done using the TCT Greeness index, Wetness index and NDVI. The pre and post were differentiated (Fig. 6) and a threshold was applied in the resulting image upon the basis that the 0 value represents pixels with no to least change (Fig.7). The results were then stacked together to produce a change detection map.
2.4. Image classification and post processing
For image classification, the 3rd PC layer together with the 24 Oct TCT wetness index and NDVI layer were all classified using the ISO (Interactive Self-Organizing Data Analysis) cluster unsupervised classification technique. It is an iterative procedure that assigns arbitrary cluster centers, classifies the pixels, calculates the cluster mean and covariance until the minimum change between iterations is achieved . The classification was done under 10 classes condition and the results were then sieved and clumped in order to remove the ’salt and pepper’ effect. The resulting images were then filtered to only retain the necessary information, such that the wetness index was for the stream, the 3rd PC was for both the boulder to pebble and finer sediments groups and the NDVI was for combined deposits (Fig. 8). These were then stacked to produce a classified map of the debris flow deposits.
An accuracy assessment was done by comparing the results of the classified map to field data that was gathered from 28-30 October. The kappa value was also calculated in order to show the extent of how the result correctly represented each classes .
3. Results and Discussion
3.1. Change Detection
Figure 9 shows the result for change detection in a portion of the study area. The decrease in NDVI value perfectly outlines the debris flow deposits since NDVI is a good tool for demonstrating a drastic change from vegetation to bare soil material. On the other hand, the decrease in TCT wetness delineated the ’passive state’ stream. This decrease was brought by the replacement of the ”wet” river by the relatively drier deposits. The decrease in TCT greeness on the other hand, indicates the areas which were previously ’greener’ but are now replaced by materials that are significantly less toned. This decrease is very prominent especially in the streams at elevated regions and portions of alluvial fans. This was probably because those areas showed the highest drastic change. The drastic change could be attributed to the fact that those areas have the highest volume of ’brown’ deposits – since they are closest to the source materials – that covered the ’green’ pixels (stream banks and rice paddies).
Some ’salt and pepper’ specs are still observable especially in the elevated regions and those that are on the fringes of clouds and cloud shadows. Though cloud masking was done, fringes and cloud shadows were still a problem. Refinement of cloud masking and detection in cloud shadows should be done in future works. The specs in the elevated regions could be easily removed by masking the elevated regions through slope selection (so as not to eliminate streams where the source materials came from). This could be done so that the results would only focus on the flatter areas where the debris flow deposits are located. However, masking was not done in this study because the primary goal is to develop an effective, semi-automated method with minimal human intervention.
Some changes weren’t also detected but this could be attributed to the threshold that was applied for each image difference. The threshold used were the most conservative thresholds based on visual inspection. However, for future work other qualitative methods for designating an unbiased threshold should be considered .
3.2. Image Classification
For image classification, only the post-typhoon image was used. Figure 10 is the result of the stacking of the unsupervised classification results for TCT wetness, 3rd PC and NDVI, as well as the field sample locations that were used for accuracy assessment. The TCT wetness index was used to delineate the stream after the typhoon. Since the satellite image used was taken six days after the event, it could be assumed that this stream delineation depicts a drier condition and is therefore a closer delineation of the streams’ ”passive state”, though if compared to the pre-typhoon image, there’s still a very obvious stream widening. Flooding on the lower right side was also detected. This could be attributed to a flooded sag pond that is the product of the presence of the Philippine Fault.
Since the 3rd PC component showed the most variation within the debris flow deposits, it was used to differentiate between boulder to pebble deposits and finer sediments. The boulder to pebble deposits are very well delineated in the alluvial fan regions. This was expected because coarser sediments brought by a debris flow event are usually deposited in the proximal region of an alluvial fan. However, there are coarse deposits that could be observed in the upper left side of the image on plain region. This result could be interpreted as a representation of the smaller coarse deposits (cobble to pebble deposits) which could be carried further away from the alluvial fan. Future efforts could include further refinement of distinguishing between larger and smaller coarse deposits.
Meanwhile, the NDVI result is very dominant in the flatter regions of the study area. They could be interpreted as combination of different deposit types including other detritus such as uprooted trees and even murky or still wet deposits. Since NDVI merely marks a drastic change from vegetation, it is safe to assume that the change was brought by a wide variety of materials. However, it is very apparent in the elevated regions that NDVI was able to detect the landslides that deposited the source material for the debris flows.
Possible refinement of the result could be done by again considering only the flatter regions though slope masking. However it was not done because of the same reason that was expressed above.
3.3. Accuracy Assessment
The accuracy assessment (Tab.3) was calculated using field data that were gathered by a team of geologists 10 days after the typhoon. The time difference between the image acquisition and the field data acquisition may have a small effect but based on the calculation of the accuracy assessment and kappa value, the results seem to be favorable and valid. The highest user accuracy was for the boulder to pebble sediments while the highest producer accuracy as achieved for both the finer sediments and the combined deposits (Tab.4). The kappa statistics shows the relevance of the method, grading it 65% better than random or chance.
Overall, the method proposed here appears to be useful for the characterization of debris flow deposits even though the method was devised in such a way that it will be semiautomatic and involves minimal intervention. NDVI was very useful for the delineation of the rapid change in the satellite image and it was therefore useful for delineating the overall debris flow deposit outline. The TCT wetness and greeness indices were very useful for change detection especially concerning the vegetated-to-flooded areas. The 3rd PC component was useful for distinguishing between debris flow deposits.
There is, however, still a problem regarding cloud shadows and cloud masking. Also, errors were also observed in the elevated regions. The problem with cloud mask and shadows could be resolved in future work by incorporating more rigorous techniques for masking and identification under shadowed areas. Errors in the elevated regions could be removed by masking through slope selection. This way, the results are limited on only the flatter areas where deposition occurred.
Future work should involve usage of more quantified thresholds as well as other types of classification methods including supervised techniques. There should also be future endeavors to apply the proposed methods on other known debris flow events to further verify its viability.
The author would like to thank the Nueva Ecija debris flow team for allowing her to use their field data and pictures.
 Philippine Atmospheric, Geophysical and Astronomical Services Administration, Severe weather bulletin #1 tropical cyclone alert: Tropical storm ”lando”, Website (October 2015). URL http://www.pagasa.dost.gov.ph/index.php/tropical- cyclones/weather-bulletin/165-tropical-cyclones/ severe-weather-bulletin/lando-2015-bulletin/1800-lando
 Joint Typhoon Warning Center, Prognostic reasoning for super typhoon 24w (koppu) warning nr 19”, Website (October 2015). URL http://weather.noaa.gov/pub/data/raw/wd/wdpn31. pgtw. txt
 National Disaster Risk Reduction and Management Council, Ndrrmc update severe weather bulletin no. 12 re typhoon ”lando”, Tech. rep., National Disaster Risk Reduction and Management Council (2015).
 Japan Meteorological Agency, Rsmc tropical cyclone advisory 171800, Website (October 2015). URL http://weather.noaa.gov/pub/data/raw/wt/wtpq20.rjtd..txt
 National Disaster Risk Reduction and Management Council, Ndrrmc update severe weather bulletin no. 13 re typhoon ”lando”, Tech. rep., National Disaster Risk Reduction and Management Council (2015).
 National Disaster Risk Reduction and Management Council, Ndrrmc update sitrep no. 22 re preparedness measures and effects of typhoon ”lando”” (i.n. koppu), Tech. rep., National Disaster Risk Reduction and Management Council (2015).
 Philippine News Agency, Pagasa plans to decommission ’lando’ from list of typhoon names, Online News Article (October 2015). URL http://www.mb.com.ph/pagasa-plans-to-decommission- lando-from-list-of-typhoon-names/
 ABS-CBN News, ’gabaldon a disaster that did not happen’, scientist says, Online news article (October 2015). URL http://rp1.abs-cbnnews.com/focus/10/22/15/gabaldon-disaster-did-not-happen-scientist-says
 D. J. Varnes, Landslide types and processes, Highway Research Board Special Report (29).
 T. C. Blair, J. G. McPherson, Alluvial fan processes and forms, in: Geomorphology of desert environments, Springer, 1994, pp. 354–402.
 R. L. Bates, J. A. Jackson, et al., Dictionary of geological terms, Vol. 584, Anchor Books, 1984.
 N. Santangelo, A. Santo, G. Di Crescenzo, G. Foscari, V. Liuzza, S. Sciarrotta, V. Scorpio, Flood susceptibility assessment in a highly urbanized alluvial fan: the case study of sala consilina (southern italy), Nat. Hazards Earth Syst. Sci 11 (2011) 2765–2780.
 J. F. Hubert, A. J. Filipov, Debris-flow deposits in alluvial fans on the west flank of the white mountains, owens valley, california, usa, Sedimentary Geology 61 (3) (1989) 177–205.
 T. Takahashi, Mechanical characteristics of debris flow, Journal of the Hydraulics Division 104 (8) (1978) 1153–1169.
 R. M. Iverson, The physics of debris flows, Reviews of geophysics 35 (3) (1997) 245–296.
 Y. K. Sohn, C. W. Rhee, B. C. Kim, Debris flow and hyperconcentrated flood-flow deposits in an alluvial fan, northwestern part of the cretaceous yongdong basin, central korea, The Journal of Geology 107 (1) (1999) 111–132.
 E. Blissenbach, Geology of alluvial fans in semiarid regions, Geological Society of America Bulletin 65 (2) (1954) 175–190.
 J. Nichol, M.Wong, Satellite remote sensing for detailed landslide inventories using change detection and image fusion, International Journal of Remote Sensing 26 (9) (2005) 1913–1926.
 D. N. Petley, W. D. Crick, A. Hart, The use of satellite imagery in landslide studies in high mountain areas, in: Proceedings of the Asian conference on remote sensing, 2002.
 G. Metternicht, L. Hurni, R. Gogu, Remote sensing of landslides: An analysis of the potential contribution to geo-spatial systems for hazard assessment in mountainous environments, Remote sensing of Environment 98 (2) (2005) 284–303.
 T. R. Martha, N. Kerle, V. Jetten, C. J. van Westen, K. V. Kumar, Characterising spectral, spatial and morphometric properties of landslides for semi-automatic detection using object-oriented methods, Geomorphology 116 (1) (2010) 24–36.
 A. Mondini, F. Guzzetti, P. Reichenbach, M. Rossi, M. Cardinali, F. Ardizzone, Semi-automatic recognition and mapping of rainfall induced shallow landslides using optical satellite images, Remote Sensing of Environment 115 (7) (2011) 1743–1757.
 K. Cheng, C. Wei, S. Chang, Locating landslides using multi-temporal satellite images, Advances in Space Research 33 (3) (2004) 296–301.
 USGS, Using the usgs landsat 8 product, Website (August 2015). URL http://landsat.usgs.gov/Landsat8_Using_Product.php
 C. Song, C. E. Woodcock, K. C. Seto, M. P. Lenney, S. A. Macomber, Classification and change detection using landsat tm data: when and how to correct atmospheric e_ects?, Remote sensing of Environment 75 (2) (2001) 230–244.
 P. S. Chavez, Image-based atmospheric corrections-revisited and improved, Photogrammetric engineering and remote sensing 62 (9) (1996) 1025–1035.
 USGS, Landsat 8 quality assessment band, Website (2014 December). URL http://landsat.usgs.gov/L8QualityAssessmentBand.php
 J. Richards, Thematic mapping from multitemporal image data using the principal components transformation, Remote Sensing of Environment 16 (1) (1984) 35–46.
 I. Jolliffe, Principal component analysis, Wiley Online Library, 2002.
 R. J. Kauth, G. Thomas, The tasselled cap–a graphic description of the spectral-temporal development of agricultural crops as seen by landsat, in: LARS Symposia, 1976, p. 159.
 D. J. Hayes, S. A. Sader, Change detection techniques for monitoring forest clearing and regrowth in a tropical moist forest, Photogrammetric Engineering and Remote Sensing 67 (9) (2001) 1067–1075.
 A. Ahmad, S. F. Sufahani, Analysis of landsat 5 tm data of malaysian land covers using isodata clustering technique, in: Applied Electromagnetics (APACE), 2012 IEEE Asia-Pacific Conference on, IEEE, 2012, pp. 92–97.
 M. L. McHugh, Interrater reliability: the kappa statistic, Biochemia Medica 22 (3) (2012) 276–282.
 P. A. Rogerson, Change detection thresholds for remotely sensed images, Journal of Geographical Systems 4 (1) (2002) 85–97.