Landslide site delineation from geometric signatures derived with the Hilbert–Huang transform for cases in Southern Taiwan

2.1 Study area

Three sites with four landslides in southern Taiwan, namely, D004, D014, D044, and D047, are selected as examples. The topographic features of the landslides include crown, main scarp, head, minor scarp, gully, reverse slope, and colluvium (Figure 1). D004, located in Meinong district, Kaohsiung City, shows a valley-like topography, covering about 0.42 km², with an average slope of 15.6° (Figure 1a). The crown has multiple ridges with several arc-shaped scarps in the middle of the landslide. The upper part of the landslide site is a flat slope, and the lower concave part of the landslide is covered by colluvial deposits, i.e., the deposit of a landslide with certain run-out distance downward from the original source of a previous landslide. The rocks at D004 are composed of Miocene sandstone, shale, and their interbeds. The colluvial deposits in the area are originated from paleo shallow landslides induced by gully erosion. In the future, deep-seated landslides with circular failure surfaces may occur in the weathered rocks or colluvial deposits as shown in Figure 4a.

Figure 1

Maps of landslide features and locations of the study areas, (a) D004, (b) D014, and (c) D044 and (d) D047. The profiles of A–A′ and B–B′ were identified by the Forestry Bureau, Taiwan in their analysis. C–C′ profile is additionally selected for this study.

Figure 2

D004 A–A′ profile and its IMFs from the original EMD. (a) The profile and its IMFs. IMF9 represents one residual component. (b) The waveforms of the first and the second IMFs. It can be seen that there are some asymmetric waveforms and two adjacent extrema which do not cross zero crossing.

D014, located in Shanlin district, Kaohsiung City, exhibits a dip-slope and valley-like topography, covering about 0.475 km², with an average slope of 25.2° (Figure 1b). The middle of the landslide deposits has several arc-shaped scarps and two significant reverse slopes. This presents a feature of gravitational slope deformation. The main scarp has an upward-sign shape. As the body of the landslide deposits is incised by a gully, the area is divided into the east and west sides [22]. The profile shape of the scarp is generally an arc. The rocks at D014 are composed of Miocene sandstone, shale, and their interbeds. Deep-seated landslides with plane failure surfaces along beddings and circular failure surfaces in the colluvial deposits may occur in the future as shown in Figure 5a.

Figure 3

D004 A–A′ elevation profile and its IMFs from the improved EMD. (a) The elevation profile and its IMFs. IMF10 represents one residual component. (b) The waveforms of the first and the second IMFs. The results show the waveforms have symmetric characteristics and that there is always zero crossing between two adjacent extrema.

D044 and D047 are located adjacently in Jiasian district, Kaohsiung City (Figure 1c). D044 covers about 0.67 km², with an average slope of 20°. The topography can be divided into three parts: the upper part is a steep ridge with a dip-slope rock, the middle is a flat slope with a scarp of jointed sandstone, and the lower part has a clear scarp with steeper terrain. D047 covers about 0.84 km², with an average slope of 17°. The topography can be divided into two parts: the upper part a dip-slope rock and the lower part a stepped-flat slope surface. The surface of the scarp is an arc shape with multiple ridges. The characteristics of the local failure are apparent, with the toe clearly delineated as a river terrace. The rocks at D044 and D047 are composed of Miocene sandstone, shale, and their interbeds. Yang et al. [23] presented the landslide as a sliding phenomenon with a buckling-induced rockslide model in dip-slope of interbedded sandstone and shale to explain the geomorphological evolutionary process of these two landslide areas. The failure of the slope in D004 was probably due to nonplanar sliding in weathered rocks (Figure 6a). The failure of the slopes in D014, D044, and D047 was likely due to a plane slide. It is possible that this plane slide failure was preceded by buckling of rock strata.

Figure 4

(a) The geomorphological features of profile D044 A–A′ annotated in black come from the profile report of the Forestry Bureau, and in red for those from the map report. (b) The IMF1 and the IMF2 spectra, respectively. (c) The spectra assembled from IMF1 and IMF2. The height is in green with an arbitrary scale.

2.2 The HHT

The implementation of HHT is based on EMD. EMD decomposes the data series into a finite number of IMFs and trend functions, with a mean value. Each of the IMF outcomes satisfies two conditions: (1) in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero [11].

To verify whether the IMF meets this basic definition, the following sifting process is taken:

1. Find the extrema of the original signal.

2. Use the cubic spline to find separately the extrema including the defined upper envelope of local maximum u_k(t) and the defined lower envelope of local minimum l_k(t).

3. Compute the local mean of envelope m_k(t) = [u_k(t) + l_k(t)]/2 to take out the component

   hk(t)=x(t)−mk(t).

4. Repeat steps 1–3 until h_k(t) can be satisfied by the definition of IMF recording c_n(t) = h_k(t).

5. Calculate the residual r_n(t) = x(t) − c_n(t).

6. If r_n(t) is identified as a trend component, the algorithm is stopped. This study takes the standard deviation computed from the two consecutive sifting results of components as the stop criterion; otherwise, steps 1–3 need to be repeated in order to find other IMFs.

From the above process, the original signal x(t) can be interpreted as the combination of n IMF components and a mean trend component r_n(t) according to equation (1)

Since this IMF is about twice the periodic characteristic of the previous one [24], a condition is included with the stop criteria of the EMD process in addition to the standard deviation: both must be satisfied. In practical experiments of this study, the threshold for standard deviation was set to 10⁻⁶, and the periodic ratio limited to a range of 1.5–2.5.

However, on the process of the elevation profile from the landslide site by HHT, we found that IMFs of the elevation data could not meet our expectation. It seemed that the existence of noises in the data might have influenced the process. Taking the IMFs of D004 AA′ elevation profile as an example (Figure 2a), we decomposed the profile into nine IMFs with the last one the residual component (IMF9 in Figure 2a). The sequence from IMF1 to IMF9 explains how the elevation is decomposed by short to long waveforms. Illustrating this problem with the first IMF (IMF1) and the second IMF (IMF2; Figure 2b), it is observed that there are some asymmetric waveforms as marked by the double-headed arrows and two adjacent extrema which do not cross the zero crossing as annotated by the single arrows.

Figure 5

(a) The geomorphological features of profile D014 A–A′ annotated in black come from the profile report of Forestry Bureau, and in red for those from the map report. (b) The IMF1 and the IMF2 spectra, respectively. (c) The spectraum assembled from IMF1 and IMF2. The height is in green with an arbitrary scale.

From the above observations, the asymmetric waveforms may be caused by the local mean in the EMD process, which is not zero. The waveforms that have not crossed the zero crossing may be affected by the extrema not being correctly estimated. In short, processing the original height series of the profile directly generates noise due to the numerical process of EMD. In order to mitigate this problem, the DI-EMD method proposed by Kopsinis and McLaughlin in 2008 [25] is introduced.

According to the sifting process step (3), the kth IMF can be obtained from the expression:

where h^(k)(t) is the temporal estimate of the kth IMF and mn(k)(t) is an estimate of the local mean of h^(k)(t) after N sifting iterations. From equation (2), it can be inferred that EMD considers the signals x(k)(t) as fast oscillations 〈h^(k)(t)〉 superimposed on slow oscillations mn(k)(t), and the sifting process aims to iteratively estimate the slow oscillating signals. As a consequence, the kth IMF is an estimate of the fast oscillating component of the signal x^(k)(t) from the signal x(t). For the estimation of the extrema, we take the first derivative of h^(k)(t), i.e., D1x(k)(t)=D1h(k)(t)+D1mn(k)(t); therefore, the extrema and the local mean of a signal can be efficiently estimated through a sifting process. That is, when a predefined number of sifting iterations is given, an estimate of D1mn(k)(t), which in turn can be subtracted from D¹x^(k)(t), is produced and leads to an estimate of the first derivative of h^(k)(t).

For the DI-EMD method, the sifting iterations used for the desired extrema estimates will be referred to as internal, and the normal EMD sifting iterations used for the current IMF estimate will be called external. Both the number of sifting iterations (it) for the desired extrema estimates and the number of sifting iterations (ex) for the original EMD are required [25]. Based on the experiments conducted in this study, the best result was achieved by two cascading DI-EMD processes. The first time, ex equals 2.5 times it; and the second time ex equals to 3.5 times it. This result is obtained in an empirical way with ratios tested ranging from a value of 1 to 3 with increments of 0.1.

Besides, to resolve the problem that two adjacent extrema (such as two adjacent local maxima or two adjacent local minima) do not have zero crossing, the EMD-based de-noising method of Kopsinis and McLaughlin in 2009 [26] is applied. The rationale is that if there are noises in the signal, when the signal is decomposed, the noises are also decomposed into finite IMFs. That is, the signal x(t) is equal to x′(t) + s(t), where x′(t) is the noiseless signal and s(t) is the noise. Assuming that the energy of all noises is captured by IMF1 (the first IMF), then the noise variance of the IMF1 E₁ may be estimated. This would mean that the thresholds of other IMFs T_i as well as the extrema of two adjacent zero crossings may likewise be estimated. All the IMF samples that correspond to zero-crossing intervals with extremum exceeding the threshold must be smoothly reduced in order for the extremum to be reduced exactly by an amount equal to the threshold. The values of β and ρ parameters in equation (4) are specifically proposed by Flandrin et al. in 2005 [27].

This can be presented as equations (3)–(6).

In which β = 0.719 and ρ = 2.01.

Here z′ij shows the threshold interval, where r_i^j indicates the jth extrema of the ith IMF, j = 1, 2,…, (M_i − 1), and M_i indicates the zero crossing number of the ith IMF.

From the above approach, the noise from the original signal may be removed. The results are shown in Figure 3a and b, where waveforms of symmetric characteristics and zero crossing between two adjacent extrema could be observed. The IMFs have been improved to more closely meet the characteristics of the narrow band. In other words, the unwanted noises in the elevation data have been removed for subsequent spectrum analysis.

Figure 6

(a) The geomorphological features of profile D044 A–A′ annotated in black come from the profile report of the Forestry Bureau, and in red for those from the map report. (b) The IMF1 and the IMF2 spectra, respectively. (c) The spectraum assembled from IMF1 and IMF2. The height is in green with an arbitrary scale.

Finally, through the Hilbert transform, the instantaneous amplitude and frequency in the IMFs are derived and a complete time–frequency distribution of energy is obtained.

For the de-noised signal X(t), the Hilbert transform can be given as equation (7):

Equation (7) is the convolution of X(t) and 1/t; P is the Cauchy principle value, so Hilbert transform can identify the local properties of X(t). The analytical signals are obtained from the conjunction of X(t) and Y(t) as equation (8):

In which

Here a(t) indicates the instantaneous amplitude of the de-noised signal X(t), θ(t) indicates the instantaneous phase angle of X(t), and the instantaneous frequency, ω(t) of X(t) as delineated by equation (9).

The de-noised signal X(t) can also be expressed as n IMFs and a residual by the abovementioned EMD method as in equation (10).

For IMFs by Hilbert transform, the same data if expanded in Fourier representation would be in the form of equation (11).

Residual component r_n(t) can be omitted.

Equation (11) can be expressed in three dimensions on the amplitude and frequency as functions of time, showing the distribution of amplitude on the frequency and time, known as the Hilbert amplitude spectrum, h(ω,t). In addition, on the time–frequency plane, it can show the amplitude or energy of the original signal to form the high or low relief of the surface.

2.3 Data preparation and processing scheme

The steps for preparing the data with HHT analysis are as follows:

1. Delimit the interested region from the 1 × 1 m digital elevation model based on the results of field geology.

2. Determine the location of the likely sliding direction of the landslide. The profile of parallel sliding direction is chosen first. The cross profile is perpendicular to the sliding direction, roughly passing the center of the landslide mass. The height profiles are generated from DEM both along and across this direction.

3. After preparing the data, HHT analysis is performed in the Matlab environment.

4. HHT and the signal de-noising analysis include three parts: (1) decompose signals at different scales to get IMFs, (2) remove the noise in elevation data and then reconstruct the de-noised elevation signals, and (3) undergo HHT to obtain the Hilbert spectrum.

The scheme of the proposed process is carried out step by step as follows:

1. collect 1 × 1 m DEM from LiDAR;

2. construct 1-m resolution discrete data points for elevation within the profile of interest;

3. decompose the discretized series of elevation data into IMFs via EMD;

4. de-noise the elevation data by reconstructing the de-noised elevation data;

5. pick the first two IMFs, i.e., IMF1 and IMF2; and

6. obtain the Hilbert spectrum of IMF1 and IMF2, individually, by the Hilbert transformation.

3.1 Elevation profile of landslide, D004 A–A′

The length of this profile is about 950 m. The landslide features from the main scarp to the toe are followed by two minor scarps, soil and weathered stratum, gully, minor scarp, colluvium, and gully (Figure 4a). The colluvium in the profile ranges horizontally from about 470 to 910 m. The distance between the main scarp and the head is about 60 m, in which there are the sag and reverse slope.

The landslide features annotated in Figure 4b are identified from the energy distribution patterns of the spectrum, where relatively obvious changes in energy are consistent with the location of the landslide features, as shown in Table 1. Generally, the Gaussian filter is used to enhance the spectrum for identifying features. In the process, the center point of each spectral feature is used for the correspondence analysis.

The spectral features observed from IMF1 and IMF2 of A–A′ (Figure 4b) with stronger energy (>0.1) are summarized as follows:

1. The energy–frequency variation from the main scarp, sag, reverse slope to the head shows a concave distribution.

2. The toe is bound by the colluvium and the gully. The energy–frequency variation near the distance from 890 to 910 m is relatively strong.

3. The upper boundary of the colluvium is about 470 m and divides the zone of depletion and the zone of accumulation. The frequency varies from 0.04 to 0.09 and the energy distribution presents a band-shaped pattern.

4. For the gully, the frequencies of 0.24 and 0.26 and the energy concentration show a sparse pattern.

Figure 4c is reassembled from IMF1 and IMF2 spectra. It appears that there are three notable zones with relatively strong energy in the spectra: the first zone is near the main scarp, the second zone is on the middle of the slope, and the third zone is on the toe. The energy between the middle of the slope and toe is higher than the energy between the main scarp and the middle of the slope, indicating the depleted mass and accumulation of the landslide body. In addition, the energy–frequency distribution from the main scarp to the head has a concave upward shape.

3.2 Elevation profile of landslide, D014 A–A′

The length of the profile is about 800 m. Figure 5a shows the landslide features of the main scarp, three minor scarps, colluvium, three gullies, and weathered stratum on the profile. In addition, from Figure 1b, it can be seen that the distance between the main scarp and the head is about 15 m, and approximately 45 m below the head of the colluvial deposit. The colluvium in the profile ranges from about 350 to 680 m bound by the gully. The minor scarp located about 580 m in the profile seems to be the active scarp (in red color in Figure 5a). The toe is bound by the terrace.

Table 2 shows the landslide characteristics corresponding to the strong energy zones (>0.05) in Figure 5b. Observation are as follows: (1) the frequency at the main scarp varies from 0.1 down to 0.052, and the frequency at the head varies from 0.09 down to 0.025 (IMF2 spectrum in Figure 5b). (2) The energies for the sag and reverse slope are relatively high in this spectrum. (3) At the toe, on the adjacent terrace, the energy is also relatively strong. (4) The frequency for the minor scarp is about 0.05. (5) For the junction of gully and colluvium at about 460 and 680 m in the spectrum along the horizontal axis, the frequency variation from about 0.05 to 0.12 and energy distribution shows a band-shaped pattern. (6) For the gully, the frequencies of 0.2 and 0.24 and the energy concentration present a parse pattern.

In Figure 5c, there are three distributed zones with relatively strong energy level, including the main scarp, the colluvium (350–680 m), and the toe. The energy between the main scarp and the reverse slope is higher than the energy between the middle and the toe. In addition, the distance between the main scarp and head is approximately 15 m, so the energy–frequency distribution presents a slight concave upward pattern.

3.3 Elevation profile of landslide, D044 A–A′

The length of the profile is about 1,150 m. Figure 6a shows the landslide features of the main scarp, two minor scarps, two colluviums, and weathered stratum on the profile. Measured from Figure 1c, the distance between the main scarp and the head is about 210 m, and the head is bound by the colluvium. The range of the two colluvia is about 210–480 m and 550–850 m, respectively. There is a gully at the toe.

The correspondence observed between the IMF1 and IMF2 spectra to landslide characteristics is as follows (Table 3 and Figure 6b): (1) the highest energy of the spectrum is located at the toe, which is adjacent to the river terrace. (2) The energy of the main scarp is greater than that of the surrounding area. The frequency at the main scarp varies from about 0.03 to 0.06. (3) The energy of the minor scarp, which is about 90 m behind the head, is very high. The phenomenon of high energy could be explained by the findings of Yang et al. [23]: this minor scarp is formed by the deposit of the run-out rock mass after plane sliding over the original slope surface. (4) In the boundary of the colluvium, which is located at a horizontal distance from 470 to 480 m and 840 to 850 m in the profile, the frequency varies from 0.025 to 0.11 and the spectrum of the energy shows a band-shaped pattern.

In Figure 6c, there are three distributed zones with relatively high energy in the energy–frequency spectra, including the main scarp, minor scarp, and the colluvium with a distance of about 600 m to the toe. Due to the distance of 210 m between the main scarp and the head, the energy distribution is not obvious, but the energy–frequency distribution between the main scarp and head also shows a concave upward pattern.

3.4 Elevation profile of landslide, D047 B–B′

This profile is about 1,450 m long. Figure 7a shows the landslide features of the main scarp, two minor scarps, sag, reverse slope, two colluviums, and weathered stratum on the profile. Figure 1c shows that the distance between the main scarp and the head is about 15 m, and there are two gullies. In addition, there are two parts of the colluvium, ranging from about 400 to 540 m and from 670 to 1,320 m in the horizontal axis of the profile.

The landslide characteristics identified and corresponding spectral features are summarized in Table 4 and presented in Figure 7b. In sum, the following observations are made: (1) the energy from the highest point of the profile goes down to the head, showing a concave upward trend in IMF1 (Figure 7b-1). But the characteristic of this concave upward trend is not significant in IMF2 (Figure 7b-2), likely resulting from the short distance between the main scarp and head, which is only 15 m. (2) At the toe, the energy is relatively strong, in contrast to the boundary of the terrace. (3) The frequency at the minor scarp is 0.06. (4) For the junction of gully and colluvium at the distance of 790 m, the frequency variation from about 0.04 to 0.1 shows a band-shaped pattern. (5) For the gully, the frequencies from 0.15 to 0.24 and the energy concentration present a sparse pattern.

In Figure 7c, it can be seen that the energy–frequency distribution presents two zones, including the main scarp to the reverse slope, and the colluvium at the distance of about 970 m to the toe. The energy distribution for the reverse slope is the strongest in the spectrum.