Towards Automated Real-Time Detection and Location of Large-Scale Landslides through Seismic Waveform Back Projection

Rainfall-triggered landslides are one of the most deadly natural hazards in many regions. Seismic recordings have been used to examine source mechanisms and to develop monitoring systems of landslides. We present a semiautomatic algorithm for detecting and locating landslide events using both broadband and short-period recordings and have successfully applied our system to landslides in Taiwan. Compared to local earthquake recordings, the recordings of landslides usually show longer durations and lack distinctive P and S wave arrivals; therefore, the back projection method is adopted for the landslide detection and location. To identify the potential landslide events, the seismic recordings are band-passed from 1 to 3Hz and the spectrum pattern in the time-frequency domain is used to distinguish landslides from other types of seismic sources based upon carefully selected empirical criteria. Satellite images before and after the detected and located landslide events are used for final confirmation. Our landslide detection and spatial-temporal location system could potentially benefit the establishment of rainfall-triggered landslide forecast models and provide more reliable constraints for physics-based landslide modeling. The accumulating seismic recordings of landslide events could be used as a training dataset for machine learning techniques, which will allow us to fully automate our system in the near future.


Introduction
Landslides and debris flows are one of serious hazards in Taiwan. Many events occurred by heavy rainfalls during rainy and typhoon seasons. The constant gravitational pull may cause downslope movement of large volumes of material. These mass movement events occur throughout the world, and under certain circumstances, some of them can produce catastrophic consequences. The general public may sometimes underestimate landslide hazards because many of the landslide events are not reported separately but are included with the descriptions of the events that have triggered them, such as volcanic eruptions, earthquakes, and major rainstorms. Globally, landslides are directly responsible for billions of dollars in damage and thousands of fatalities and injuries each year (e.g., [1]).
Situated at the boundary between the Philippine Sea Plate to the east and the Eurasian Plate to the west, Taiwan is a mountainous and seismically active region (e.g., [2]). About one-third of Taiwan's land area is covered by rugged mountainous topography. Temperature differences between the East Asian continent and the Pacific Ocean drive monsoonal flow that carries warm, moist air from the Indian Ocean and the Pacific Ocean and delivers periodic monsoon rains and several typhoons each year. The combination of rugged mountainous topography, abundant water, and frequent earthquakes makes Taiwan particularly susceptible to landslide hazards. The 2009 Typhoon Morakot produced an accumulated rainfall of 2777 mm and triggered over 400 large-scale landslides involving over 100,000 m 2 area throughout southern Taiwan, leaving nearly 700 fatalities and over $4.7 billion USD in damage (e.g., [3]). It is therefore crucial to establish the capability to detect and locate hazardous landslides in (near) real-time in Taiwan for emergency response and hazard mitigation purposes.
For rainfall-induced landslides, surveillance photographs from aircrafts and/or optical remote-sensing satellites are not immediately available often due to extreme weather conditions. Previous studies have shown that the landslide can be represented as an equivalent shallow seismic source with either a nearly horizontal single-force mechanism or a horizontal thrust fault mechanism (e.g., [4][5][6][7]). These results suggest that it is possible to adapt seismic techniques for (near) real-time earthquake detection and epicenter location to detect and locate landslides rapidly.
However, such an adaptation is not straightforward. Complications arise from two major differences between the seismograms generated by a landslide and those generated by a typical tectonic earthquake. First, seismograms from landslides often lack clear arrivals of P and S body wave phases, which are key signatures for identifying and locating earthquakes in the majority of (near) real-time earthquake detection/location techniques. Second, seismograms generated by landslides are often highly complex and it is not yet clear what may constitute a robust waveform signature for distinguishing landslides from other seismic events.
It appears that landslides are capable of generating sufficient low-frequency energies that can be detected by broadband seismic stations located at large distances. Inversions of such long-period waveforms (0.02~0.05 Hz) allowed Lin et al. [8] to recover the locations and source mechanisms of 52 landslides and submarine slumps during the 2009 Typhoon Morakot. More recently, the near real-time "landquake" monitoring system of Chao et al. [9] is also based upon the inversion of low-frequency waveforms (0.02~0.05 Hz) in their detection module. However, automated detection and inversion algorithms based purely upon low-frequency waveforms have the risk of misidentifying medium-sized and large earthquakes, which can also generate significant low-frequency energies, such as landslides, and may miss some landslides that do not generate sufficient low-frequency energies.
Our semiautomated landslide detection and location system mainly relies upon seismic waveform data in the frequency range 1~3 Hz. We use low-frequency waveforms (e.g., in the 0.02~0.05 Hz band) mainly for the purpose of rejecting teleseismic events, which may have similar timefrequency spectrum as landslides in the 1~3 Hz band. This choice of the frequency band allows us to utilize waveforms not only from the broadband seismic network but also from the extensive short-period seismic network in Taiwan. The increased station density reduces station gaps and improves azimuthal coverage, which is often a crucial factor in reducing location uncertainties.
Because of the lack of clear P and S arrivals on landslide seismograms, in our automated location algorithm, we have adopted the back projection algorithm, originally developed for locating nonvolcanic tremors and glacial earthquakes, for locating landslides in Taiwan (e.g., [10,11]). The travel times used in back projection are computed using 3D ray tracing in one of the latest 3D seismic velocity models of Taiwan. To improve the efficiency of our landslide location code, we compute all travel time data in advance and archive them on disk. The code for locating the landslide only needs to read the data from the disk without performing timeconsuming ray-tracing calculations.
We have tested our automated system on past continuous seismic recordings from broadband and short-period seismic networks in Taiwan. Our semiautomated system correctly identified the majority of the landslide events and provided the locations with~2 km uncertainty in less than 5 minutes of the event origin time. When applied to real-time continuous seismic recordings, our semiautomated system may help reduce the emergency response time in disastrous landslide events.
Furthermore, to understand the complex mechanisms of rainfall-triggered landslides, the knowledge from different disciplines, such as geology, hydrology, geomorphology, and soil sciences, need to be considered in the system. The initiation time and location of a landslide event are critical for understanding the threshold conditions of landslide occurrences. The landslide events determined based on seismic recordings provide accurate initiation times and reliable locations, and this type of information could benefit other rainfall-triggered landslide studies. For example, accurate rainfall-triggered landslide catalog and rainfall records could be very helpful in developing a more reliable rainfallbased landslide forecast model for hazard management and mitigation (e.g., [12][13][14]). In addition, accurate occurrence times of rainfall-triggered landslides could provide more reliable conditions (e.g., rainfall duration, rainfall intensity, and cumulative rainfall) for physics-based landslide modeling (e.g., [15,16]).

Data and Landslide
Waveform Characteristics 2.1. Seismic Networks in Taiwan. In Taiwan, the majority of previous seismic studies of landslides (e.g., [1][2][3][4]) have been using only waveform recordings from the Broadband Array in Taiwan for Seismology (BATS), which is operated by the Institute of Earth Sciences (IES) of Academia Sinica and the Central Weather Bureau (CWB) [2]. In addition to BATS, the CWB also operates a telemetered, digital, short-period seismic network ( Figure 1) composed of 1 Hz Teledyne S13 seismometers [17]. For the frequency band used in this study (i.e., 1~3 Hz), waveforms recorded by the short-period network can also be utilized to improve landslide detection accuracy and to constrain landslide locations. By including the stations from the short-period network, in addition to the stations in the BATS, into our landslide detection/ location system, we reduced the average interstation distance from over~30 km to less than~20 km ( Figure 1). A subset of the short-period stations is colocated with the broadband stations in the BATS. In Figure 2, we show the high waveform similarities within the frequency band 2 Geofluids 1~3 Hz of a typical landslide recorded by colocated broadband and short-period stations. Our back projection algorithm uses the envelopes of the 1~3 Hz waveforms. In Figure 2(f), we show that the envelope of the short-period recording within the 1~3 Hz frequency band is practically identical to that of the broadband recording in the same frequency band with a correlation coefficient larger thañ 0.9. In our study, waveform recordings from both networks are utilized and, at those locations where broadband and short-period instruments are colocated, the waveform data from the broadband instruments take precedence.

Waveform
Signature of the Landslide in the Time-Frequency Domain. Previous studies have shown that within the 1~3 Hz frequency band, waveforms of the landslide have a distinctive time-frequency domain signature that can be used to distinguish them from local tectonic earthquakes and teleseismic earthquakes (e.g., [9,[18][19][20][21][22][23]). Figure 3 shows seismic waveform recordings of a landslide ( Figure 3, left), a typical local earthquake (Figure 3, center), and a typical teleseismic event (Figure 3 right) recorded by Taiwan's seismic network. Compared to the local earthquake signal, the landslide signal lacks clear onsets of seismic body wave arrivals (i.e., P and S wave arrivals) and has a longer temporal duration. On the time-frequency spectrum, the landslide event shows a distinctive isosceles triangle or shield pattern (Figure 3 bottom left), which is substantially different from a right-angle triangle pattern ( Figure 3, bottom center) typical for local earthquakes. The comparison of time-frequency domain spectrums suggests that high-frequency energies generated by landslides arrive later and decay earlier within the duration of the signal than those generated by a typical local tectonic earthquake. This observation is consistent with results from previous studies, and we also use this observation to distinguish landslides from local earthquakes in our landslide detection algorithm.
However, the time-frequency domain spectrum patterns are not sufficient to distinguish teleseismic events from landslides. Certain portions of the teleseismic recording ( Figure 3, bottom right) can show a similar time-frequency domain spectrum pattern within the 1~3 Hz frequency band as that of a typical landslide (Figure 3 bottom center). In our landslide detection algorithm, we use low-frequency waveforms (0.02~0.05 Hz) to distinguish teleseismic events from landslides. At such low frequencies, teleseismic recordings at different stations tend to have high waveform similarities and the differential travel time at neighboring stations tends to be very small. These low-frequency features of teleseismic events are mainly due to the high incidence angles of the incoming wave and can be used as additional criteria for rejecting teleseismic earthquakes.

Methods and Representative Examples
Conceptually, we can organize our landslide detection/ location workflow sequentially into 5 different stages ( Currently, these 5 stages are carried out iteratively in a semiautomated process. In the following, we will discuss each stage separately.

Seismic Data Preprocessing.
The main purpose of the data preprocessing stage is to enhance the waveform signature of landslides while suppressing contributions from local tectonic earthquakes and other unwanted seismic events. This stage has 5 separate steps, and all the steps were designed to facilitate the detection and location algorithms downstream in our flowchart ( Figure 4). First, we cut all continuous raw waveform recordings from both the broadband and the short-period seismic networks into 70-minute-long segments with 10-minute-long overlaps ( Figure 5(a)). The choices of time segment length and overlap length are configurable and can be adjusted manually to optimize the overall performance of the entire system. Second, we remove the mean and the linear trend of each time segment and then apply a minimal-phase Butterworth bandpass filter with corner frequencies at 1 and 3 Hz ( Figure 5(b)).  Figure 1: Distributions of broadband (gray symbols) and shortperiod (red squares) stations in Taiwan. The background color shows regional topography and bathymetry.

Geofluids
Third, we calculate the root mean square (RMS) of threecomponent recordings ( Figure 5(c)). In general, the landslide signal has time durations longer than about 15 seconds. To suppress short-term fluctuations, we smooth the RMS by computing its moving average with a 10-second-long time window ( Figure 5(d)).
Fourth, the smoothed RMS is normalized using a threshold value chosen specifically for each individual time segment. An example is shown in Figure 5(d), in which we can see clearly an event with the largest amplitude located at 12 minutes, which is due to a local M L 3.6 earthquake, and an event with a smaller amplitude located at~57 minutes, which is due to a landslide. In this example, the threshold value is set at the 99th percentile of the entire time segment. Values above this threshold are clipped, and values below this threshold are normalized using the threshold value. Figure 5(d) shows the smoothed RMS after the threshold normalization (red line) in comparison with the one before the threshold normalization (black line). The event associated with the landslide at~57 minutes is enhanced, while the event associated with the local earthquake at~12 minutes is suppressed. We have tried different threshold values for the normalization, and generally values between the 99th and 99.5th percentiles of the smoothed RMS amplitudes gave reasonable results for our landslide detection and location stages later on.
Finally, we compute the envelope of the time segment obtained from the fourth step through Hilbert transform and feed the result to the detection algorithm in the next stage.
3.2. Landslide Event Detection. The majority of conventional earthquake-location techniques mainly rely upon picking the onset time of body wave arrivals, such as P and S waves (e.g., [24]). The seismic recordings of landslides usually lack clear onset time of body wave phases ( 4 Geofluids (e.g., [19]). The back projection method for locating seismic sources does not rely upon precise picks of body wave onset times and has been successfully applied to locate tremors, as well as many different types of seismic sources (e.g., [19,[25][26][27]). In this study, we have adopted the back projection method to detect and locate landslides in Taiwan. We back project a segment of the envelope function obtained at the end of stage 1 (section 3.1) and perform a weighted stack across multiple stations according to the following equation: where E i is the normalized envelop at station i; t i is the model-predicted travel-time from a source located at x with origin time t to station i; τ 0 , τ 1 is the time window containing the normalized and shifted envelopes used for stacking; W i is the weighting factor for station i; and N is the total number of stations used in the back projection. The name "back projection" comes from the temporal shift t i x, t , which propagates the envelop function backwards in time from the station location to the space-time coordinate x, t of the potential source. If the seismic source is in fact located at x and occurred at time t, the normalized and shifted envelopes at different stations E i will be systemically aligned and stack constructively, which would lead to a local maximum of the stacking function S at the source space-time coordinate x, t ; otherwise, they will stack incoherently and will not produce a maximum. By detecting the emergence of such maxima, we can detect the occurrences and estimate the locations of landslides and other seismic events.
In the landslide detection stage, we use a spatial grid with 3 km grid spacing that covers the surface topography of Taiwan at elevations higher than 10 m above the sea level. Previous studies have suggested that the dominant seismic energies generated by landslides are primarily Rayleigh waves and/or shear waves (e.g., [9,[18][19][20][21][22][23]). We therefore calculated and archived the S wave travel times from each spatial grid point to all stations in BATS and the short-period network using one of the latest 3D seismic velocity models of Taiwan [28]. Currently we use the 3D ray tracing software, simulPS12 [29], for travel time calculations. Surface waves usually arrive slightly later than S waves; we therefore multiply the model-predicted S wave travel times by 1.08 to obtain estimates of surface wave travel times for the shifting operation.
Seismic signals generated by landslides usually have longer duration than local tectonic earthquakes. Increasing the length of the time window τ 0 , τ 1 used in the stacking operation can potentially suppress contributions from unwanted local tectonic earthquakes. We have tried various time window lengths, and in general, a length of about 20 seconds provides reasonable estimates of landslide sources.
For (near) real-time applications, the processing speed of the detection stage is critically important. To accelerate the process of detection and location and to reduce the computational cost, currently we only use the 10~20 At the detection stage, we assign the same weighting factor W i to the envelope function of each selected station. The local maxima in the stacking function S indicate potential landslide events (Figures 5(e) and 5(f)). In this study, we select the local maxima with values larger than 6 times the median absolute deviation (MAD) of the stacking function as the space-time coordinates of potential landslide events. The detection criteria could be adjusted based on the researchers' actual applications.
3.3. Landslide Event Inspection. In general, due to their short temporal durations, we can effectively reject the majority of small local tectonic earthquakes (M L < 2.5) during stage 2. The primary purpose of stage 3 is to reject any local earthquakes that we failed to reject during stage 2 and any teleseismic events with substantial energy in the 1~3 Hz frequency band.
For every event detected in stage 2, we cut the seismic waveform recordings from 30 s before the origin time to 180 s after the origin time at all available broadband and short-period stations. The selected waveform data are then converted into the time-frequency domain through continuous wavelet transform, and we check the patterns of the time-frequency spectrum. As shown in Figure 3, seismograms generated by landslides usually show an isosceles triangle or shield spectrum pattern in the time-frequency domain, while the spectrums of local earthquakes often show a right-angle triangle pattern. We use this criterion to reject any remaining local earthquakes from the set of events detected during stage 2.
Teleseismic events are rejected based upon a different criterion, since their time-frequency spectrum may have similar patterns as those of landslides in the 1~3 Hz frequency band. Seismograms from teleseismic events usually contain significant low-frequency energies, and the waveforms have very high similarities among the stations. We apply a Butterworth bandpass filter with corner frequencies at 0.02 and 0.05 Hz to the broadband recordings. For waveforms with a signal-to-noise ratio larger than 5, we then calculate the waveform correlation coefficients among station pairs and also the differential delay times. If the average of the correlation coefficients is higher than a given threshold (e.g., 0.8 in our current implementation)  6 Geofluids and the differential delay times between station pairs are much smaller than the expected travel times for the S wave, we reject the event.
An example is shown in Figure 6. During the 2015 Typhoon Soudelor, one event detected during stage 2 has landslide-like waveforms at the high-frequency band (1~3 Hz), but it was rejected in stage 3 based upon our teleseismic rejection criteria. The waveforms at the lowfrequency band (0.02~0.05 Hz) show very high similarities ( Figure 6(b)). The differential delay times between station pairs are much smaller than the corresponding travel times of S arrivals. For example, the model-predicted S travel time differences from the detected source location are 39.88 s for MASB and SSLB, 26.32 s for SSLB and YHNB, and 26.99 s for MASB and YULB. But the differential delay times for the above station pairs are only 3.75 s, 1.25 s, and 0.75 s, respectively. The significant variations in time differences among the station pairs indicate that the lowfrequency energies were generated by a teleseismic event, which was verified by the 1000-second long low-frequency recordings (Figure 6(c)).

Landslide Event Relocation and Magnitude Estimates.
All the events that have survived stage 3 are then relocated at stage 4. At this relocation stage, the seismic recordings are selected based upon a more stringent signal-to-noise ratio (SNR) criterion and the source location search is performed on a refined grid surrounding the preliminary location obtained from stage 2. Figure 7 shows examples of the relocation stage for landslide events during the 2012 Typhoon Saola and the 2015 Typhoon Soudelor periods. The inputs at this stage are the envelopes of the 3-component RMS of 1~3 Hz bandpass-filtered broadband and short-period recordings (Figures 7(a) and 7(b), right). Different from stage 2, each envelope is normalized by the maximum amplitude at the station. For each envelope, we select a time window from 5 s before the model-predicted S arrival time to 25 s after At this stage, we search for the optimal source location on a 60 km wide square grid with 1 km grid spacing and centered at the location obtained from stage 2 (Figures 7(a) and 7(b), left). We use the SNR of each envelope as the weighting factor, W i , and the length of stacking time window τ 0 , τ 1 is 5 s long. The maximum of the stacking function over the refined grid gives the estimated landslide location of stage 4 (Figures 7(a) and 7(b), left).
The amplitudes of recorded ground motions could be used to quantify the energy released by sources. Similar concepts have been applied to landslide magnitude estimate in previous studies [19,30]. For example, Lin et al. [30] use the long-period signals (20-50 s) of large landslides occurred during 2009 Typhoon Morakot to define an empirical equation for landslide magnitudes. Kao et al. [19] used the maximum amplitudes from the recordings for magnitude estimates. In this study, we adopted the definition of local magnitude for earthquakes in Taiwan area [31] to estimate the magnitude of a detected landslide. In the local magnitude formula, the distances from the event are based on our refined landslide location to individual stations. The magnitudes of the landslides which occurred during 2012 Typhoon Saola and 2015 Typhoon Soudelor are M L = 1 6 and M L = 2 0, separately.

Landslide Event Confirmation.
To confirm the landslides detected and located during stages 1-4, we manually check the satellite images before and after each typhoon event. If a significant landslide could be identified within 10 km of our estimated landslide location, we consider our detection and location system to be successful and then catalog the landslide event. If our system did not detect the landslide event 9 Geofluids or our estimated location is more than 10 km away from the actual location, we consider this case as a failure. Figure 8 shows the satellite images of two landslide events detected and located by our system (Figure 7) during the 2012 Typhoon Saola and the 2015 Typhoon Soudelor periods. The 2012 landslide event is located by our system at 121.51°E, 24.62°N. The location of this landslide identified from the satellite images is around 121.57°E, 24.57°N. The location difference is about 8.6 km. For the 2015 landslide event, the optimal location estimated using our system is at 120.81°E, 23.21°N and the closest event identified from the satellite image is at around 120.78°E, 23.21°N. The location difference is about 3.3 km. Since the SNRs of the seismic waveforms generated by the 2012 landslide event are relatively low (Figure 7(a)), as compared with those from the 2015 event (Figure 7(b)), it is expected that the location error for the 2012 event is larger than that of the 2015 event.

Comparison with Landslide Events Determined by Other
Methods. We have applied our system to detect and locate landslides during several previous typhoon events, including the 2009 Typhoon Morakot. In previous studies [18,19], only the broadband seismograms recorded by the BATS were used for landslide locations. We have included all available broadband and short-period seismograms from CWB's networks in our workflow. The landslide events listed in [18,19] have all been successfully detected and located by our system. In addition to those landslide events, our system also identified a few more landslide events that were missed by previous studies. Figure 9 shows some additional landslide examples detected by our algorithm and confirmed by satellite images. In general, locations determined by our system are consistent with those obtained from satellite images.

Other (Near) Real-Time Automated Landslide Detection/ Location Systems in Taiwan.
At the current stage, the majority of previous landslide studies in Taiwan have not yet considered real-time automation in detail (e.g., [8,18,19]). The latest (near) real-time automated landslide detection/ location system proposed in [9] has a detection module that is based primarily upon low-frequency (0.02~0.05 Hz) seismic signals recorded by broadband instruments only. The landslide events detected by this system have been cataloged and archived at the website http://collab.cv.nctu.edu.tw/ catalog.html. A close examination of this catalog shows that, without human intervention, a considerable number of local earthquakes have been misclassified as landslides (i.e., false positive). Taiwan is a region with high seismicity, and the frequent medium-sized local earthquakes can generate considerable low-frequency energies. If the local earthquake has a shallow hypocenter and a low-angle thrust mechanism, there is a high probability of misidentifying such an earthquake as a landslide.
A second drawback of using only low-frequency signals for landslide detection is that there is a high possibility of missing landslides that do not generate sufficient lowfrequency energies (i.e., false negative). In Figure 10, we show broadband recordings of two landslides, one during the 2012 Typhoon Saola (Figure 9, left) and the other during the 2015 Typhoon Soudelor (Figure 9, right). For both events, the waveforms at the 1~3 Hz frequency bands show significant energies and good SNR, but at the 0.02~0.05 Hz band, the Locations of these two landslides determined using our seismic system are shown in Figure 7. Geofluids 2012 landslide shows poor SNR. Our semiautomated landslide detection/location system is based primarily upon the broadband and short-period recordings in the 1~3 Hz brand.

Landslide Relocation Algorithms.
In [9], the landslide relocation algorithm is based upon the cross-correlation of the horizontal envelope function at neighboring stations. Such a cross-correlation-based relocation algorithm has been successfully applied to relocate nonvolcanic tremors [32], whose seismograms usually lack clear body wave onsets just like the seismograms generated by landslides. However, we must point out that there are substantial differences between the source mechanisms of nonvolcanic tremors and those of typical landslides. For nonvolcanic tremors, weak sources with similar mechanisms occur in a short period of time (e.g. a few minutes) and the time intervals between significant energy bursts (e.g., S waves) at neighboring stations are almost identical. A direct consequence of such source mechanisms is that the envelope functions at neighboring stations have very high similarities (e.g., Figure 7 in [33]); therefore, the cross-correlation lag times between stations can be treated as the differential S wave travel times for relocating tremors. The source mechanisms of landslides can be very different from those of nonvolcanic tremors, and the similarities of the envelope functions at neighboring stations are not as high, especially when taking into account the relatively larger station gaps of the broadband network in Taiwan. Examples of highly dissimilar envelope functions from a landslide are shown in Figure 7. Those fundamental differences in source mechanisms may increase the errors in the landslide locations obtained from the cross-correlation relocation algorithm.

Conclusions and Future Developments
We have developed a semiautomatic system for landslide detection and location using broadband and short-period seismic recordings, and we have successfully applied our system to landslides in Taiwan. The design of our algorithm takes into account the possibility that some landslides 11 Geofluids may not generate significant low-frequency energies, so we mainly use seismic recordings in the high-frequency band (1~3 Hz) for detection and location. To prevent earthquakes (local or teleseismic) from being classified as landslides, we have implemented rejection criteria based upon time-frequency domain spectrum patterns and lowfrequency (0.02~0.05 Hz) interstation waveform correlations. The satellite images before and after the detected and located landslide events near the affected regions are used for landslide event confirmation.
Mass-movement hazards, such as landslides, are one of the most dangerous natural disasters in the world (e.g., [34]). Many countries, including Taiwan, suffer significant fatalities and damages from rainfall-triggered landslides each year during monsoon/typhoon seasons. Accurate estimates of occurrence times and locations of landslide events are playing an increasingly important role in improving our understanding of landslide initiation conditions, physicsbased modeling, and developing forecast models. We expect our landslide detection/location system to make valuable contributions to the global efforts in landslide hazard analysis and mitigation.
Recent advances in deep neural networks (e.g., [35]) have opened up opportunities to fully automate our landslide detection/location system. Successes in image classification using deep convolutional neural networks (CNN) suggest that we can adopt CNN for distinguishing seismograms generated by landslides from those generated by earthquakes using time-frequency domain spectrograms automatically.
The grid search for finding optimal source locations and origin times can be accelerated, perhaps significantly, by adopting massively parallel coprocessors such as the graphic processing unit (GPU). The use of the GPU for accelerating the detection of small-to medium-sized earthquakes through cross-correlation-based template matching has been successfully implemented and tested [36]. Similar GPU acceleration techniques can also be adopted for landslide detection and location. We will leave those studies for future publications.

Data Availability
The waveform data used for this study are from the Geophysical Database Management System (GDMS) in Taiwan operated by the Central Weather Bureau (CWB) and the Broadband Array in Taiwan for Seismology (BATS) operated by the Institute of Earth Sciences, Academia Sinica (IES). The satellite images are provided by the Center for Space and Remote Sensing Research at the National Central University.