Detection of Earthquake-Induced Landslides during the 2018 Kumamoto Earthquake Using Multitemporal Airborne Lidar Data

: A series of earthquakes hit Kumamoto Prefecture, Japan, continuously over a period of two days in April 2016. The earthquakes caused many landslides and numerous surface ruptures. In this study, two sets of the pre- and post-event airborne Lidar data were applied to detect landslides along the Futagawa fault. First, the horizontal displacements caused by the crustal displacements were removed by a subpixel registration. Then, the vertical displacements were calculated by averaging the vertical di ﬀ erences in 100-m grids. The erosions and depositions in the corrected vertical di ﬀ erences were extracted using the thresholding method. Slope information was applied to remove the vertical di ﬀ erences caused by collapsed buildings. Then, the linked depositions were identiﬁed from the erosions according to the aspect information. Finally, the erosion and its linked deposition were identiﬁed as a landslide. The results were veriﬁed using truth data from ﬁeld surveys and image interpretation. Both the pair of digital surface models acquired over a short period and the pair of digital terrain models acquired over a 10-year period showed good potential for detecting 70% of landslides.


Introduction
An Mw 6.2 earthquake hit the Kumamoto Prefecture in Kyushu Island, Japan, on 14 April 2016, at 21:26 (JST). The epicenter was located in the Hinagu fault with a shallow depth. Twenty-eight hours later, another Mw 7.0 earthquake occurred in the Futagawa fault on April 16 at 01:25 (JST), close to the Hinagu fault. The first event was called a "foreshock" and the second one the "main shock". Figure 1a shows the locations of these causative faults and the Japanese national GNSS Earth Observation Network System (GEONET) stations in the source area [1]. 75-cm horizontal movements to the east-northeast were observed at the Kumamoto station, while 97-cm horizontal movements to the southwest were recorded at the Choyo station during the main shock. Strong shaking, level 7 on the Japan Meteorological Agency (JMA) seismic intensity scale was observed in Mashiki Town in both the foreshock and the main shock events. A total of fifty-five people died in the earthquake sequence, mostly due to the collapse of wooden houses in Mashiki town and landslides in Minami-Aso village. 8673 buildings collapsed, and approximately 35,000 buildings were severely damaged [2]. The earthquake also caused the interruption of lifelines, e.g., gas, electric power and water. Road and railway networks were affected by landslides and surface ruptures, making it difficult to access the affected area soon after the earthquake.
provide a different kind of information about landslides [30,31]. Zhang and Lu [32] reviewed the application of remote sensing to landslides, and Jaboyedoff et al. [33] summarized the use of airborne Lidar and ground Lidar in landslide investigations. A high-resolution DEM obtained by airborne Lidar measurement can detect different types of landslides; however, it is uncommon due to its high cost. The 2016 Kumamoto Earthquake was a rare case when three temporal Lidar flights were carried out within one month. Thus, detecting landslides using those precious data sets is a valuable research. In this study, the detection of landslides that occurred in Mashiki town and Nishihara village, Kumamoto Prefecture, Japan, were conducted using two sets of Lidar data acquired before and after the 2016 Kumamoto Earthquake. An effective method of landslide detection is proposed. The Remote sensing technologies have been employed to understand the extent and the degree of various damages after natural disasters [3][4][5][6][7][8][9][10][11]. For the 2016 Kumamoto Earthquake, many studies have been conducted to assess earthquake-induced damage using remote sensing data [12][13][14][15]. The authors calculated the pre-and co-event coherence from three temporal ALOS-2 PALSAR-2 images to classify building damage levels in Mashiki town [16]. One day after the foreshock (15 April 2016), Asia Air Survey Co., Ltd. carried out a Lidar surveying flight from 15:00-17:00 (JST), to record the effects of the earthquake. One week after the main shock, they flew a second mission from 10:00-12:00 (JST) on 23 April. The authors' group used those two temporal Lidar data to estimate the coseismic crustal displacements and extract debris from the collapsed buildings in Mashiki town [17,18]. Three weeks after the main shock, Geospatial Information Authority of Japan (GSI) also conducted another Lidar surveying flight on 8 May 2016. The third flight covered the largest area from the east of Mt. Aso to the Minami ward of Kumamoto City. The authors attempted to extract collapsed buildings for Mashiki town from this Lidar data and another Lidar data acquired 10 years ago [19].
Remote sensing techniques for mapping landslides can be classified into two types. One method is to extract surface changes using optical images or synthetic aperture radar (SAR) intensity images [20][21][22]. Optical images can provide details of landslides, but they are affected by sun and weather conditions. SAR intensity images can overcome this problem, but their results are influenced by slope orientation. SAR intensity images are usually combined with optical images to extract landslides [23]. The other method is to measure surface movements using InSAR analysis or Lidar data [24][25][26][27]. InSAR analysis can detect landslides by two different processes. The differential InSAR method or stacking InSAR methods (PSI, SBAS) are good at monitoring slow-moving landslides [28,29]. For rapid debris flows, the significant surface changes indicate the decorrelation in the InSAR pairs. A digital elevation model (DEM) generated from InSAR analysis or measured by airborne Lidar can provide a different kind of information about landslides [30,31]. Zhang and Lu [32] reviewed the application of remote sensing to landslides, and Jaboyedoff et al. [33] summarized the use of airborne Lidar and ground Lidar in landslide investigations. A high-resolution DEM obtained by airborne Lidar measurement can detect different types of landslides; however, it is uncommon due to its high cost. The 2016 Kumamoto Earthquake was a rare case when three temporal Lidar flights were carried out within one month. Thus, detecting landslides using those precious data sets is a valuable research.
In this study, the detection of landslides that occurred in Mashiki town and Nishihara village, Kumamoto Prefecture, Japan, were conducted using two sets of Lidar data acquired before and after the 2016 Kumamoto Earthquake. An effective method of landslide detection is proposed. The obtained results were verified by comparing with the results of visual interpretation. The objective of this study is not only to detect landslides but also to compare the performance using two different sets of Lidar data.

Landslides Associated to the 2016 Kumamoto Earthquake in the Study Area
In the main shock, upper-6 JMA seismic intensity was observed in the Minami-Aso village. The strong shaking caused many landslides surrounding the Aso volcano and in areas within 10-km from the Futagawa fault. The National Research Institute for Earth Science and Disaster Resilience (NIED) released several versions of a coseismic landslide inventory [34]. In the latest edition on 27 June 2016, 1744 landslides were traced from the 460-km 2 interpretation range. GSI identified the location of 82 landslides larger than 0.01 km 2 and 661 landslides larger than 0.005 km 2 by visual interpretation using the aerial photos taken on 16, 19 and 20 April 2016 [1]. The number of landslides caused by rainfall and aftershocks was updated to 1044, using the aerial photos taken in July 2016.
During our field surveys, many landslides were observed. Two ground photos are shown in Figure 1b. One is the largest landslide that occurred in the Minami-Aso village. This landslide was 200 m in width and 700 m in length. The maximum depth of erosion reached 20 m. The landslide overwhelmed the Bungo highway and destroyed the Aso-Ohashi Bridge. This photo was taken by one of the authors on 9 August 2016. In addition to those at the Aso volcano, many small landslides occurred along the Futagawa fault. The other ground photo shows a landslide in Dozono district, Mashiki town. This photo was taken on 7 June 2016. The trees covering the hill slid down with soils along the steep slope.
Many studies of landslides have been conducted on from the 2018 Kumamoto Earthquake. Doi et al. [35] used geological and geophysical methods to understand the generation mechanism of the landslides that occurred in Minami-Aso village. Tajima et al. [36] conducted landform change analysis of the central cones of the Aso volcano using satellite and aerial images. Tamkuan and Nagi [37] combined Landsat-8 and interferometric ALOS-2 coherence data to classify building damage, liquefaction and landslides. Kim and Kim [38] used object-based classification with the truth color RapidEye images to extract the landslides in Minami-Aso village. Uemoto et al. [39] identified the outliers of the largest landslide based on height and amplitude differences using airborne X-band SAR data. Most of those studies focused on landslides surrounding the Aso volcano. The authors attempted to detect landslides close to the Futagawa fault using the Lidar data acquired on 15 and 23 April 2016 [40,41].
Our study area is located at the boundary of Mashiki town and Nishihara village, along the Futagawa fault. The study area is shown in Figure 2 by three red frames. Their identified IDs are 2KD672, 2KD673 and 2KD762. Each frame covered an area of 2 km × 1.5 km. According to the previous study [17], more than 2-m right-lateral strike-slip and 2-m subsides were detected on the west side of the fault in this area. Due to these large crustal movements, many small landslides and numerous surface ruptures occurred. were provided as orthorectified truth color images with 16-cm spacing, whereas the photos on 23 April were provided as orthorectified truth color images with 23-cm spacing.

Figure 2.
Enlarged area of the white frame in Figure 1a. The coverage of the DSMs acquired by Asia Air Survey Co., Ltd. is enclosed by the dashed line, and the coverage of the DSMs acquired by GSI is shown in the black gridlines. The vertical differences between the pre-and post-event DSMs acquired by Asia Air Survey Co., Ltd., after removing the influence of the crustal displacements, also overlaps the base map. Our target areas are shown with red frames.
The other Lidar pair were acquired on 15 May 2006, and 8 May 2016, by GSI. Since the purpose of the pre-event Lidar data was to create a 5-m digital terrain model (DTM), the point density was low (0.72 points/m 2 ). The post-event Lidar was taken to detect damage, with a high density of 5.93 points/m 2 . The coverage area is shown in Figure 2 with black gridlines. Both the original random point data including ground and surface objects, and the processed DTMs with 2 m spacing were provided. The original random point data were used to detect the collapsed buildings in the center of Mashiki town [19]. Aerial photos were also acquired during the flights. The aerial photos obtained in 2006 were panchromatic images with 30-cm spacing, whereas the photos acquired in 2016 were truth color images with 50-cm spacing.
The timeline of the 2016 Kumamoto earthquake event and the acquisitions of the Lidar data is shown in Figure 3. All the Lidar data were transformed to raster images by ENVI Lidar software. Considering the different point densities, the data acquired by Asia Air Survey Co., Ltd. were resampled with 1-m spacing, whereas the data acquired by GSI were resampled with 2-m spacing (pre-event) and 1-m spacing (post-event). Examples of a pre-event aerial photo, the DSM acquired on    Figure 1a. The coverage of the DSMs acquired by Asia Air Survey Co., Ltd. is enclosed by the dashed line, and the coverage of the DSMs acquired by GSI is shown in the black gridlines. The vertical differences between the pre-and post-event DSMs acquired by Asia Air Survey Co., Ltd., after removing the influence of the crustal displacements, also overlaps the base map. Our target areas are shown with red frames.

Lidar Data and Preprocessing Analysis
Two sets of pre-and post-event airborne Lidar data were used in this study. One pair was the Lidar data acquired on 15 April and 23 April 2016, by Asia Air Survey Co., Ltd. [42]. The average point density on 15 April was 1.5-2 points/m 2 , and that on 23 April was 3-4 points/m 2 . Both sets of Lidar data were acquired using a Leica ALS50II instrument and the same pilot and airplane. The raw point clouds have been resampled in 50-cm meshes. The common coverage areas are shown in Figure 2 by the dashed line. This data set had been used to estimate the coseismic displacements and detect collapsed buildings [17,18]. The vertical differences between the pre-and post-event data after removing the influence of the crustal displacements is also shown in Figure 2. In our target area along the Futagawa fault, several significant differences can be observed, which are considered to be landslides. The authors attempted to detect landslides in the target area using a simple thresholding value (2 m) [40,41]. However, the results included much noise and did not evaluate well. In the flights, aerial photos were acquired at the same locations as the Lidar data. The aerial photos on 15 April were provided as orthorectified truth color images with 16-cm spacing, whereas the photos on 23 April were provided as orthorectified truth color images with 23-cm spacing.
The other Lidar pair were acquired on 15 May 2006, and 8 May 2016, by GSI. Since the purpose of the pre-event Lidar data was to create a 5-m digital terrain model (DTM), the point density was low (0.72 points/m 2 ). The post-event Lidar was taken to detect damage, with a high density of 5.93 points/m 2 . The coverage area is shown in Figure 2 with black gridlines. Both the original random point data including ground and surface objects, and the processed DTMs with 2 m spacing were provided. The original random point data were used to detect the collapsed buildings in the center of Mashiki town [19]. Aerial photos were also acquired during the flights. The aerial photos obtained in 2006 were panchromatic images with 30-cm spacing, whereas the photos acquired in 2016 were truth color images with 50-cm spacing.
The timeline of the 2016 Kumamoto earthquake event and the acquisitions of the Lidar data is shown in Figure 3. All the Lidar data were transformed to raster images by ENVI Lidar software. Considering the different point densities, the data acquired by Asia Air Survey Co., Ltd. were resampled with 1-m spacing, whereas the data acquired by GSI were resampled with 2-m spacing (pre-event) and in 2006 were panchromatic images with 30-cm spacing, whereas the photos acquired in 2016 were truth color images with 50-cm spacing.
The timeline of the 2016 Kumamoto earthquake event and the acquisitions of the Lidar data is shown in Figure 3. All the Lidar data were transformed to raster images by ENVI Lidar software. Considering the different point densities, the data acquired by Asia Air Survey Co., Ltd. were resampled with 1-m spacing, whereas the data acquired by GSI were resampled with 2-m spacing (pre-event) and 1-m spacing (post-event). Examples of a pre-event aerial photo, the DSM acquired on 15 April 2016, and the DSM and the DTM acquired on 15 May 2006, are shown in Figure 4, respectively. The figures show the images in grid 2KD673. The elevation of this area ranged from 80 m to 250 m. The southeast side had higher altitude than the northwest side. Comparing the DSMs in 2016 (b) and 2006 (c), several land-use changes could be confirmed during the 10-year period. In the DTM (d), buildings and trees had been removed, and the fine relief could be clearly observed.   Since the main shock caused more than 2-m horizontal displacements and 2-m subsidence in the study area [17], a subpixel-based registration was adopted to remove the horizontal movements. The pre-event images were set as the base. Tie points were generated automatically using the crosscorrelation method. The maximum allowable error was set as less than 0.5 pixel. To improve the accuracy of the DTM registration, tie points obtained using the pre-and post-event DSMs were applied to the DTMs. After the registration, the horizontal differences were less than 40 cm. All the post-event images were warped in 1-m spacing. Then, the vertical displacements were estimated by averaging the vertical differences in the 10 × 10 pixel (10 m × 10 m) windows. The obtained vertical displacements using the DSMs and the DTMs are shown in Figure   Since the main shock caused more than 2-m horizontal displacements and 2-m subsidence in the study area [17], a subpixel-based registration was adopted to remove the horizontal movements. The pre-event images were set as the base. Tie points were generated automatically using the cross-correlation method. The maximum allowable error was set as less than 0.5 pixel. To improve the accuracy of the DTM registration, tie points obtained using the pre-and post-event DSMs were applied to the DTMs. After the registration, the horizontal differences were less than 40 cm. All the post-event images were warped in 1-m spacing. Then, the vertical displacements were estimated by averaging the vertical differences in the 10 × 10 pixel (10 m × 10 m) windows. The obtained vertical displacements using the DSMs and the DTMs are shown in Figure 5 Finally, the corrected vertical differences after removing both the horizontal and vertical coesismic displacements were obtained from the DSMs acquired by Asia Air Survey Co., Ltd. and the DTMs acquired by GSI. They are shown in Figure 6. Several negative regions (blue colors) can be identified along the fault from a, whereas only one significant changed region could be confirmed from b. Many red lines (increased height) were observed parallel to the fault. These lines fit the boundaries of trees at the different altitudes. The height increase may be caused by the different Lidar point locations in the two flights. The blue lines are the locations of felled trees, which were confirmed from the 23 April 2016, aerial photo. In the southeast part of Figure 6b, several changed areas can be observed far from the fault. Comparing the aerial photos taken in 2006 and 2016, those regions were confirmed to be artificial land-cover changes.

Methodology
As mentioned in Section 2, we attempted to detect landslides using the DSMs acquired on 15 and 23 April 2016 [40,41]. Using a simple thresholding value of 2 m, the landslides were detected with much noise. Thus, we improved the method by adding the topographic features. The flowchart of the proposed method is shown in Figure 7. First, the significantly increased and decreased regions were extracted by the thresholding method. Then, the maximum slope angle in each decreased region Finally, the corrected vertical differences after removing both the horizontal and vertical coesismic displacements were obtained from the DSMs acquired by Asia Air Survey Co., Ltd. and the DTMs acquired by GSI. They are shown in Figure 6

Methodology
As mentioned in Section 2, we attempted to detect landslides using the DSMs acquired on 15 and 23 April 2016 [40,41]. Using a simple thresholding value of 2 m, the landslides were detected with much noise. Thus, we improved the method by adding the topographic features. The flowchart of the proposed method is shown in Figure 7. First, the significantly increased and decreased regions were extracted by the thresholding method. Then, the maximum slope angle in each decreased region

Methodology
As mentioned in Section 2, we attempted to detect landslides using the DSMs acquired on 15 and 23 April 2016 [40,41]. Using a simple thresholding value of 2 m, the landslides were detected with much noise. Thus, we improved the method by adding the topographic features. The flowchart of the proposed method is shown in Figure 7. First, the significantly increased and decreased regions were extracted by the thresholding method. Then, the maximum slope angle in each decreased region (erosion) was calculated. The regions of decrease on flat ground (collapsed buildings or moving vehicles) were removed. The increased regions (depositions) close to the decrease regions on the slope were identified according to the slope aspect. Finally, a hole-filling step and a closing filter were applied to clean up the results.

Thresholding Method
Thresholding is the simplest way to extract a region of interest. However, determining the threshold value is the most important and difficult step. There are many methods to find an optimal threshold value [43][44][45]. To simplify the approach and maintain the accuracy, we determined the threshold value manually by comparison with aerial photos. The average of the DSM differences between 15 and 23 April was 2 cm. The standard deviation was 2 m. In the previous study [40,41], ±2 m was used as the threshold value for the DSMs acquired by Asia Air Survey Co., Ltd. As shown in the ground photo of Figure 1b, the landslides occurred on slopes covered by high trees. The felled trees will cause a height reduction larger than 2 m. Thus, the previous results included much noise. An enlarged area within the black frame in Figure 6 is shown in Figure 8. The extracted regions by the different threshold values are overlapping on the ortho-image on 23 April 2016 (Figure 8a,b). There were two large landslides in this area, which were also identified by both NIED and GSI. Those two landslides could be completely detected from the DSMs, even using the strictest threshold value (−6 m). The threshold value of −3 m extracted several tree changes. Thus, we used −4 m as the best threshold value. In the DTM results, the shapes of the two landslides could not be extracted completely when the threshold value was set as −4 m. Furthermore, the threshold value of −2 m extracted many non-landslide regions. The threshold value was determined to be −3 m for the DTMs.

Thresholding Method
Thresholding is the simplest way to extract a region of interest. However, determining the threshold value is the most important and difficult step. There are many methods to find an optimal threshold value [43][44][45]. To simplify the approach and maintain the accuracy, we determined the threshold value manually by comparison with aerial photos. The average of the DSM differences between 15 and 23 April was 2 cm. The standard deviation was 2 m. In the previous study [40,41], ±2 m was used as the threshold value for the DSMs acquired by Asia Air Survey Co., Ltd. As shown in the ground photo of Figure 1b, the landslides occurred on slopes covered by high trees. The felled trees will cause a height reduction larger than 2 m. Thus, the previous results included much noise. In this study, we tried the thresholding values of An enlarged area within the black frame in Figure 6 is shown in Figure 8. The extracted regions by the different threshold values are overlapping on the ortho-image on 23 April 2016 (Figure 8a,b). There were two large landslides in this area, which were also identified by both NIED and GSI. Those two landslides could be completely detected from the DSMs, even using the strictest threshold value (−6 m). The threshold value of −3 m extracted several tree changes. Thus, we used −4 m as the best threshold value. In the DTM results, the shapes of the two landslides could not be extracted ground changes were difficult to observe as they were obstructed by trees. Since there are no trees in the DTMs, these two small erosions have a high probability of being landslides. Similar to the results from the DSMs, many depositions were extracted.

Modification Using the Topographic Features
To remove the regions of decrease caused by collapsed buildings or moving vehicles, the slope information was used to discriminate landslides. We calculated the slope angles from the DTMs in 2006 and 2016, respectively. Considering the point density of the Lidar data acquired in 2006, a 15 × 15 pixel (15 m × 15 m) window was used in the calculation. Comparing those two results, some slopes became flat after landslides. Thus, only the pre-event slope information was effective for detection. According to a previous study [47], a landslides' frequency peak falls between 15° and 40°. The high spatial resolution of the DTMs captured the fine undulations on the boundary of roads and agricultural fields. We investigated the slope angle of the boundary of those flat regions, which were up to 20°. Thus, the threshold value for dividing the sloped and the flat areas was defined as 20°. We calculated the maximum value of the slope angle within each erosion. If the maximum value was less than 20°, it would be removed from the results. One erosion extracted from the DSMs was removed in this step due to its small slope value. It is shown in Figure 8c by the white line. The extracted results from the DTMs were also modified using the same slope information.
As shown in Figure 8c,d, many depositions were extracted. We distinguished the depositions due to the landslides by linking them to the erosions. In this step, the aspect angle was used. The aspect angle was also calculated from the DTM in 2006 using a 15 × 15 pixel window. The aspect angle The extracted erosions and depositions using the optimal threshold values are shown in Figure 8c,d. From the DSMs, four large erosions were extracted, including two landslides and two felled tree regions. Three small erosions were confirmed as felled trees, whereas two were not shown in the aerial photo. Many depositions were extracted, including two landslides and other changes of tree heights. From the DTMs, four erosions, including two landslides and two other regions were obtained. The two small extracted erosions were close to the locations of the erosions extracted from the DSMs. From the aerial photo, the felled trees could be confirmed at this location; however, the ground changes were difficult to observe as they were obstructed by trees. Since there are no trees in the DTMs, these two small erosions have a high probability of being landslides. Similar to the results from the DSMs, many depositions were extracted.

Modification Using the Topographic Features
To remove the regions of decrease caused by collapsed buildings or moving vehicles, the slope information was used to discriminate landslides. We calculated the slope angles from the DTMs in 2006 and 2016, respectively. Considering the point density of the Lidar data acquired in 2006, a 15 × 15 pixel (15 m × 15 m) window was used in the calculation. Comparing those two results, some slopes became flat after landslides. Thus, only the pre-event slope information was effective for detection. According to a previous study [47], a landslides' frequency peak falls between 15 • and 40 • . The high spatial resolution of the DTMs captured the fine undulations on the boundary of roads and agricultural fields. We investigated the slope angle of the boundary of those flat regions, which were up to 20 • . Thus, the threshold value for dividing the sloped and the flat areas was defined as 20 • . We calculated the maximum value of the slope angle within each erosion. If the maximum value was less than 20 • , it would be removed from the results. One erosion extracted from the DSMs was removed in this step due to its small slope value. It is shown in Figure 8c by the white line. The extracted results from the DTMs were also modified using the same slope information.
As shown in Figure 8c,d, many depositions were extracted. We distinguished the depositions due to the landslides by linking them to the erosions. In this step, the aspect angle was used. The aspect angle was also calculated from the DTM in 2006 using a 15 × 15 pixel window. The aspect angle was 0 • at the north direction, increasing clockwise. The average value of the aspect angle for each erosion was estimated. When the average aspect angle was smaller than 45 • or larger than 315 • , the erosions shifted north to search for the associated depositions. When the angle was between 45 • and 135 • , the erosions shifted east. When the angle was between 135 • and 225 • , the erosions shifted south. When the angle was between 225 • and 315 • , the erosions shifted west. The shift was 1 m for the results from the DSMs, and 2 m for the results from the DTMs, the same as the original spacing of the raster images. If the shifted erosion overlapped on a deposition, that deposition was linked to the erosion. Only the linked depositions were counted as parts of landslides. Then, the erosions and depositions were merged into one object as a landslide. The hole in the extracted object was filled. In addition, a closing filter using a 3 × 3 pixel window was applied to clean up the shapes of the results.
The modified results are shown in Figure 9. From the DSMs, five landslides were extracted. The extracted regions of two landslides were increased by merging the depositions. Four landslides were extracted from the DTMs. The traced landslide regions by NIED are also shown in Figure 9 by the white lines. For the two recognizable landslides, the shapes matched the regions traced by NIED. The confusion matrixes were calculated in this area, using the polygon traced by NIED as the truth data. They are shown in Table 1. The producer accuracy of the DSM results was 87% and the user accuracy was 32%. Since many felled trees were extracted wrongly as landslides, the user accuracy was low. The producer accuracy of the DTM results was 77%, and the user accuracy was 56%. Compared to the DSM results, the commission errors reduced significantly. In both the results, the high producer accuracies showed the potential of the proposed method for detecting landslides. The landslides dominated only 10% of the area, thus, both the overall accuracies were over 90%. The kappa coefficient was 0.44 for the DSM results and 0.63 for the DTM results. was 32%. Since many felled trees were extracted wrongly as landslides, the user accuracy was low. The producer accuracy of the DTM results was 77%, and the user accuracy was 56%. Compared to the DSM results, the commission errors reduced significantly. In both the results, the high producer accuracies showed the potential of the proposed method for detecting landslides. The landslides dominated only 10% of the area, thus, both the overall accuracies were over 90%. The kappa coefficient was 0.44 for the DSM results and 0.63 for the DTM results.

Results and Verification
The proposed approach and the defined threshold values were applied to the entire target area in grids 2KD 672, 2KD673 and 2KD762. The obtained results are shown in Figure 10. From the DSMs acquired on April 15 and 23, 2016, 110 regions were detected to be landslides. Most of the extracted results were located close to the Futagawa fault. The total area of the detected landslides was 0.097 km 2 , including 0.066-km 2 erosions and 0.031-km 2 depositions. The volume of the results was also calculated. The decreased volume was 8.1 × 10 5 -m 3 and the increased volume was 1.7 × 10 5 -m 3 . The increased volume was 20% of the decreased one. The DSMs measure the heights at the top of trees. When a landslide occurred, the trees on the slope fell and slid with soils to the ground. The decreased volume included the space of the standing tress, whereas the increase volume included the felled trees. Thus, the decreased volume was much larger than the increased volume. A concept image of this phenomenon is shown in Figure 11  The obtained results were compared with the landslides traced by NIED. NIED investigated 21 landslides in the target area, covering a 4.5 × 10 4 m 2 area. Sixteen landslides were identified successfully from the pre-and post-event DSMs, and 12 landslides were identified from the DTMs. Two new omission errors in the DSM results were two small landslides with lower height changes. The nine omissions in the DTM results showed the limited capability on small landslide detection. The problem is mainly due to the low point density of the pre-event data. A comparison of our results and the landslides identified by NIED is summarized in Table 3. Two enlarged comparisons are shown in Figure 12b,c. Five landslides were detected from both the DSMs and DTMs, one landslide was detected only from the DSMs, and one could not be detected by either. The two landslides that could not be detected from the DTMs were smaller than other detected landslides. Comparing their shapes, the regions of piled-up soils were difficult to distinguish from the DSMs. The threshold value of the deposition was set to 2 m for the DSMs, a little high to extract accumulated soils. However, the DSMs were sensitive to the erosions of a small landslide. The threshold value of the deposition from the DSMs was set to 1 m. Thus, most of the piled-up soils were detected. The detected depositions were larger than the visual interpretation carried out by NIED. The confusion matrixes calculated in these two enlarged areas (0.4 km 2 ) are shown in Table 4. The DSMs could detect 53% of the affected regions, whereas the DTMs could detect 67%. The kappa coefficient was 0.56 for the results from the DSMs and 0.55 for the results from the DTMs, showing moderate agreement with the visual interpretation result.
According to those comparison, the two temporal DSMs acquired in a short period are sensitive for detecting landslides. However, the results were affected by trees. Felled trees were difficult to distinguish from landslides without additional information. When the two temporal Lidar data sets were taken over a long period, DTMs are more suitable for landslide detection. In this study, the low point density of the pre-event Lidar caused the low resolution of the DTM. Small landslides could not be detected by our approach. However, the proposed method would work better using highdensity Lidar data. Since the Lidar data acquired by Asia Air Survey Co., Ltd. have been resampled into grids, the DTM could not be generated. The influence of time lag on the processing of the DTMs could not revealed. The total area of the detected landslides was 0.064 km 2 , including 0.027-km 2 of erosions and 0.037-km 2 of depositions. The detected depositions were larger than the erosions. However, when they were transformed into volume, the erosion volume was 1.5 × 10 5 m 3 and the deposition volume was 1.6 × 10 5 m 3 . The total changed volume was much smaller than the results obtained from DSMs. Since the DTMs represented the ground altitude, only soil changes were detected. Thus, the decreased volume was close to the increased volume.
The locations of landslides investigated by GSI using the aerial photos taken in April 2016, are also shown in Figure 10. In the target area, 13 small landslides (larger than 0.005 km 2 ) and one large landslide (larger than 0.01 km 2 ) were indicated. Comparing with the DSM results, 11 out of 14 landslides were successfully detected. The three omitted locations were landslides that occurred after the foreshock, which was confirmed from the aerial photo. An example of two omitted locations is shown in Figure 12a. In the aerial photo taken on 15 April 2016, these two landslides had occurred. This means our results using the pre-and post-main-shock DSMs detected all the landslides caused by the main shock. Comparing with the DTM results, 10 out of 14 landslides were successfully detected. Three of the omissions showed only height increases on the post-event DTMs. Without a neighboring erosion, the depositions were not identified as landslides. Thus, they could not be detected in our results. One omission showed that the height decreased; however, the area that decreased more than 3-m was less than 100 m 2 . Those errors may be caused by the low point density of the pre-event Lidar data. The producer accuracy and user accuracy of our results counted by number is summarized in Table 2. The low user accuracy of the DSM results was mainly caused by the felled trees due to surface ruptures or soil movements, which can be partially observed in Figure 9a    The obtained results were compared with the landslides traced by NIED. NIED investigated 21 landslides in the target area, covering a 4.5 × 10 4 m 2 area. Sixteen landslides were identified successfully from the pre-and post-event DSMs, and 12 landslides were identified from the DTMs. Two new omission errors in the DSM results were two small landslides with lower height changes. The nine omissions in the DTM results showed the limited capability on small landslide detection. The problem is mainly due to the low point density of the pre-event data. A comparison of our results and the landslides identified by NIED is summarized in Table 3. Two enlarged comparisons are shown in Figure 12b,c. Five landslides were detected from both the DSMs and DTMs, one landslide was detected only from the DSMs, and one could not be detected by either. The two landslides that could not be detected from the DTMs were smaller than other detected landslides. Comparing their shapes, the regions of piled-up soils were difficult to distinguish from the DSMs. The threshold value of the deposition was set to 2 m for the DSMs, a little high to extract accumulated soils. However, the DSMs were sensitive to the erosions of a small landslide. The threshold value of the deposition from the DSMs was set to 1 m. Thus, most of the piled-up soils were detected. The detected depositions were larger than the visual interpretation carried out by NIED. The confusion matrixes calculated in these two enlarged areas (0.4 km 2 ) are shown in Table 4. The DSMs could detect 53% of the affected regions, whereas the DTMs could detect 67%. The kappa coefficient was 0.56 for the results from the DSMs and 0.55 for the results from the DTMs, showing moderate agreement with the visual interpretation result. According to those comparison, the two temporal DSMs acquired in a short period are sensitive for detecting landslides. However, the results were affected by trees. Felled trees were difficult to distinguish from landslides without additional information. When the two temporal Lidar data sets were taken over a long period, DTMs are more suitable for landslide detection. In this study, the low point density of the pre-event Lidar caused the low resolution of the DTM. Small landslides could not be detected by our approach. However, the proposed method would work better using high-density Lidar data. Since the Lidar data acquired by Asia Air Survey Co., Ltd. have been resampled into grids, the DTM could not be generated. The influence of time lag on the processing of the DTMs could not revealed.

Conclusions
The landslides caused by the 2016 Kumamoto Earthquake in Mashiki town and Nishihara village were detected using two sets of airborne Lidar data. In estimating the coseismic displacements, the DSMs acquired over a 10-year period were found to be unsuitable for detection due to the presence of significant surface changes. Thus, detection was conducted using the DSMs on 15 and 23 April 2016, and the DTMs on 15 May 2006, and 8 May 2016. From comparisons using NIED and GSI image interpretation, the locations of 71% of the landslides could be identified. The DSMs were more sensitive than the DTMs for detecting small landslides. However, the DTM results were more stable compared to the DSM results. When the time lag between the two Lidar data acquisitions is long, DTMs can obtain more promising results than DSMs. The point density is an important consideration for the accuracy of the results. Higher resolution DTMs can improve the detection accuracy. Since GSI conducts Lidar surveys on a regular basis, the proposed method can be applied when post-event Lidar data are available.