Modeling Accumulated Volume of Landslides Using Remote Sensing and DTM Data

Abstract

Landslides, like other natural hazards, such as avalanches, floods, and debris flows, may result in a lot of property damage and human casualties. The volume of landslide deposits is a key parameter for landslide studies and disaster relief. Using remote sensing and digital terrain model (DTM) data, this paper analyzes errors that can occur in calculating landslide volumes using conventional models. To improve existing models, the mechanisms and laws governing the material deposited by landslides are studied and then the mass balance principle and mass balance line are defined. Based on these ideas, a novel and improved model (Mass Balance Model, MBM) is proposed. By using a parameter called the “height adaptor”, MBM translates the volume calculation into an automatic search for the mass balance line within the scope of the landslide. Due to the use of mass balance constraints and the height adaptor, MBM is much more effective and reliable. A test of MBM was carried out for the case of a typical landslide, triggered by the Wenchuan Earthquake of 12 May 2008.

1. Introduction

Landslides are one of the most dangerous natural hazards. They often occur suddenly, destroy critical infrastructure and ecological systems, and cause large-scale damage to buildings, as well as cause human casualties and economic losses. It was reported that landslide events account for 42% [1] of the global incidences of natural disasters, with average annual economic losses of millions to billions of US dollars in such countries as Canada, Nepal, and Sweden. In China, over 300,000 sites, where landslides and other geological hazards occurred, have been investigated and it was found that landslides are becoming more frequent and serious. In the areas affected by the Wenchuan Earthquake of 12 May 2008, more than 20,000 people were killed by landslides [2]. In February 2010, massive landslides on Portugal’s Atlantic island of Madeira killed 38 people [3].

To counter the effects of landslides, engineering measures, such as retaining walls, anti-sliding piles, and anchor cables have been widely used. A key design parameter for these engineering structures is the volume of the landslide. Accurate calculation of this parameter is a challenging task, especially for those landslides with large-scale and unclear slide surfaces. In addition, the landslide volume is an important index for determining the scale, for stability analysis and risk assessment and also for evaluating the investment needed for dealing with landslides [4–8]. The landslide volume can also be used to predict secondary hazards, such as debris flows and dammed lakes. As for the morphological characteristics of a landslide, its location, distribution, surface, and area are easily acquired, whereas its volume is relatively difficult to estimate as this requires detailed geometric information concerning the complex slope failure surface, irregular sub-surfaces, and the deposit thickness and height [9,10].

There are a number of different methods that can be used to calculate the landslide volume. In some studies, the landslide volumes were progressively divided into regular or nearly regular units and were calculated by means of integral approximation approaches, for which some key parameters were acquired through field survey, geophysical prospecting, or drilling [11–13]. In another group of studies, statistical relations between landslide volume and characteristics including frequency, surface area, and the inversion of surface displacement data were established in the case of regional landslides [8,14–16]. The statistical relations between landslide volume and precipitation, especially antecedent or triggering precipitation, were also set up [17,18]. Other statistical relations between landslide volume and seismic action, including earthquake magnitude or intensity, etc., were adopted in the case of landslides triggered by earthquakes [19–21]. All the above methods are dependent on numerous samples or statistical data, which are difficult to acquire and also difficult to apply in other locations.

Thanks to the rapid development of remote sensing technologies in recent decades, remote sensing data have been widely used for landslide studies. The use of stereo photogrammetry and airborne LiDAR (Light Detection And Ranging) for the study of landslides is undergoing fast growth. Both of these techniques have been used to produce high-resolution Digital Terrain Models (DTM) or Digital Elevation Models (DEM) and can be used to detect and analyze the spatial distribution of landslides through the sliding activity, e.g., cracks, scarps, and folds [22–26]. Recent studies, based on DTMs, or combined DTM and remote sensing data of unstable terrain, focused on landslides boundary changes and on the estimated volume of accumulated or removed material [21,27–34]. In general, such combined methods calculate landslide volumes by calculating elevation differences between pre- and post-landslide DTMs.

In this paper, errors that may reduce the accuracy of commonly used model are, first, analyzed in detail and an innovative model, based on the mass balance principle, is proposed. According to this principle, the accumulated volume of the rotational and translational landslides can be acquired by automatically searching for a mass balance line. The paper also describes the theoretical basis of this work before the model is applied to the case of the 12 May 2008, Wenchuan Earthquake. The advantages and disadvantages of the proposed model are also discussed.

2. Method and Model

2.1. Problem Analysis

The model commonly used for landslide volume calculation, which compares the elevation differences between pre- and post-landslide DTMs within the landslide-affected area as identified from remote sensing images, is [35]:

$$V={\int }_a^b\mathit{Sdh}$$

(1a)

$$V=\Delta S\Delta h_1+\dots +\Delta S\Delta h_n=\Delta S\sum _{i=1}^n\Delta h_i$$

(1b)

where V (m³) is the landslide’s accumulated volume, S (m²) is the horizontal area between a and b, and h (m) is the elevation difference between the pre- and post-landslide DTMs, with a and b representing the lowest and the highest elevation respectively. Equation (1b) is a discrete form of Equation (1a), where n is the number of segments; ΔS (m²) and Δh (m) are the horizontal segment area and the elevation difference of the segment, respectively.

According to the above model (referred to as Height Difference Model, HDM), DTM subtractions allow the variations in the topography to be quantified: negative values for the Z-coordinates correspond to subsidence or ablation of rock and soil, these values can be used to derive the removed volume (see Figure 1a). Positive values reflect movement where subsidence is combined with the advance of the landslide [36] and these values can be used to derive the accumulated volume of the landslide. This model is quite simple and is relatively accurate as long reliable data are available. It has been used by some researchers [27,35,37–40] and has also been adopted by several commercial software packages, such as the MicroStation [35]. However, errors in the plane position (X- and Y-coordinates), height value (Z-coordinate), and height baseline difference may reduce the accuracy of the results produced by this model.

Errors in the pre- and post-landslide plane positions produce co-registration errors between the two corresponding DTMs. It is well known that a small movement in the plane position may correspond to a large elevation change in mountainous areas where there are many sudden changes in elevation (see the steep surface lines in Figure 1a). Consequently, this inconformity will appear in the final result, which means that the above model is quite sensitive to co-registration errors between the pre- and post-landslide DTMs. The co-registration error was noticed and addressed by using matching 3D or 2.5D data when using multi-temporal data to determine the deformation field of landslides [41–45].

The error in the DTM height values relates to many factors, such as the original data, the model, the interpolation algorithm, the mesh grid distance and so on. This error is random and always appears directly in the value of Δh in the model described above. Moreover, there is also another kind of error which directly reduces the accuracy of the value of Δh to a large degree. This error arises from the systematic differences in the height datum, the geographical reference ellipsoid, and the projection method. Therefore, both random errors and systematic errors are inevitable and produce an error in Δh in Equation (1). To perform calculations using different DTMs that had inhomogeneous baselines, [27] selected five sub-areas in order to correlate the average altitudes (The model introduced by [27] is referred to as the Advanced Height Difference Model, AHDM). The results showed that it is possible to correlate and position all of the homologous control points within areas of overlap between the DTMs, which is a crucial fact for further volume determinations. Unfortunately, apart from AHDM, there are few published articles on other improvements to the conventional model.

If the above errors did not exist then, outside the landslide-affected area, Δh would be zero and S would have no impact on V. In reality, however, these errors listed above are inevitable and the error in S will contribute to the error in V.

2.2. Modeling

The above error analysis gives a theoretical understanding of the conventional model and illustrates where there is room for improvement. Figure 1b is a typical photograph of a landslide [46]. It clearly shows the topographical changes caused by the landslide; there is also material movement, with upper material sliding or rolling down, leaving a distinct eroded area at the top. The material that fallen down piled up into a deposit. Figure 1a shows a side projection of this kind of landslide and can considered to be a simplified model. In Figure 1a, the black line stands for the post-landslide topographical profile and the grey line represents the pre-landslide profile. The difference between the two profiles corresponds to the removed volume for the part under the pre-landslide surface and to the accumulated volume for the part above the pre-landslide surface. The dotted line represents the part of the deposited volume, which is under the pre-landslide profile (“filling up” volume). In some cases, where points C and D coincide with point O, there is no filling up volume. The total deposited volume is then equal to the accumulated volume and the initial volume is equal to the removed volume.

As the pre- and post-landslide topographical profiles are obtained from the pre- and post-landslide DTMs, the commonly used model (Equation (1)) can only calculate the accumulated volume. Therefore, from here onwards, the landslide volume mentioned in this paper is always the accumulated volume that can be calculated using Equation (1).

A closer look at Figure 1 may naturally lead to the following conclusion. The landslide only moves the material rather than adding to or removing it. The accumulated material is composed of soil and stones that fallen from the upper side of the slope, which means that the amount of accumulated material should be the same as that removed from the upper slope. This equivalence, which provides a perfect chance to improve the accuracy of Equation (1), is the theoretical basis of the current study and is referred to as the mass balance principle. The relevant definitions are given below.

2.2.1. Mass Balance Principle

Landslides, by their nature, are physical movements. A landslide triggers only the movement of material—with a possible change in density—but no changes in its mass or composition. The mass of the landslide-triggered accumulated material is the same as that of the material removed from the upper slope.

According to the definition of the mass balance principle, a new concept called the mass balance line is defined as follows.

2.2.2. Mass Balance Line

Within the landslide range, the boundary along which the mass of accumulated material equals that removed from the upper slope is defined as the mass balance line.

The mass balance line is, by definition, a theoretical line. In Figure 1a, for example, point O is the projection of the balance line on to the profile. The volume above the balance line—the removed volume—equals the volume under it—the accumulated volume. According to this definition, the total deposited volume is also equal to the initial volume.

To sum up, the proposed model can be described as follows.

If m⁻ and m⁺ represent the masses of the removed material and the accumulated material respectively, then:

$$m^{-}=m^{+}$$

(2)

and as volume (V) times density (ρ) equals the mass then:

$$m^{-}=V^{-}{\rho }^{-};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{m}^{+}=V^{+}{\rho }^{+}$$

(3)

It is well known that soil usually becomes loose after movement [47]. In other words, the density of the accumulated material is lower than that of the removed material. Given the change in density, a quantity called the density adaptor, α (less than 1), is introduced here:

$$\alpha {\rho }^{-}={\rho }^{+}$$

(4)

From Equations (2–4) we can infer that:

$$V^{-}=\alpha V^{+}$$

(5)

Equation (5) contains both the volume of the accumulated material (V⁺) and that of the removed material (V⁻). According to Equation (1), both V⁻ and V⁺ are calculated from the elevation difference between the pre- and post-landslide DTM measurements. So, if the elevation difference changes, both V⁻ and V⁺ also change, and eventually V⁻ and V⁺ reach equilibrium in Equation (5). Accordingly, an improved landslide volume calculation model based on the mass balance is given as the following pseudo-code (referred to as Mass Balance Model, MBM):

$$\begin{array}{l}[1]:\hspace{0.17em}\hspace{0.17em}\Delta H_i^k=H_i^{\mathit{after}}-H_i^{\mathit{before}}+H_k^b;\\ [2]:\hspace{0.17em}\hspace{0.17em}{V}_k^{-}=\Delta S\sum _{i=1}^n\mathit{abs}\left(\Delta H_i^k\right),\Delta H_i^k<0;\\ [3]:\hspace{0.17em}\hspace{0.17em}{V}_k^{+}=\Delta S\sum _{i=1}^n\Delta H_i^k,\Delta H_i^k>0;\\ [4]:\hspace{0.17em}\hspace{0.17em}\mathit{if}\left(V_k^{-}=\alpha V_k^{+}\right)\\ [5]:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{print}\hspace{0.17em}\hspace{0.17em}{H}_k^b,V_k^{+},V_k^{-};\\ [6]:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{break};\\ [7]:\hspace{0.17em}\hspace{0.17em}\mathit{else}\\ [8]:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{H}_{k+1}^b=H_k^b+\mathit{Hstep};\\ [9]:\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{go}\hspace{0.17em}\hspace{0.17em}[1];\\ [10]:\hspace{0.17em}\hspace{0.17em}\mathit{end}\hspace{0.17em}\mathit{if}\end{array}$$

(6)

where, $H_i^{\mathit{after}}$ represents the elevation in the post-landslide DTM while $H_i^{\mathit{after}}$ represents the elevation in the pre-landslide DTM. ΔS is the area of each DTM grid. The parameter α is the density adaptor. i = 1, 2, …, n, where n is the number of discrete units in the landslide; k = 0, 1, …, m, where m is the number of loops in the calculation.

H^b, a height adaptor, is defined to adjust the pre-landslide and post-landslide DTM height data. It has two roles in MBM: first it is used to eliminate the systematic difference between the DTM elevations before and after the landslide and, secondly, in MBM, it is the driver that makes an automatic search for the balance line. In this sense, the volume calculation is translated into an automatic searching process for the balance line. This is shown in Figure 2a, which is a simple slice through a landslide. The red line A₀B₀ represents the pre-landslide surface profile, whereas the blue solid line, AB, represents the post-landslide slope profile. According to mass balance principle, triangle A₀OA represents the volume of the removed area while triangle B₀OB represents the volume of the accumulated volume. When MBM is used, if the DTM profile line is A₁B₁, and its elevation is lower than that of line AB, V⁻ is bigger than V⁺. In order to achieve the ideal mass balance, A₁B₁ needs to be moved upwards and H^b will be a positive number. As H^b increases, the gap between V⁻ and V⁺ diminishes until it reaches zero and A₁B₁ is adjusted to be exactly the same as AB. Conversely, if the DTM profile line is A₂B₂ then H^b will be negative and the adjustment will finally cause A₂B₂ to coincide with AB. The vertical distance between AB and A₁B₁ or the distance between AB and A₂B₂ is H^b. In this way, using MBM the automatic calculation of the landslide volume can be realized with a reliability that is obviously superior to that obtained using Equation (1), which has no binding conditions at all.

From the above analysis, it can be seen the range of H^b is deemed as being (−∞,+∞) with an original value of 0. The only parameter that it is required to set is the step size for H^b. This will be discussed in Section 3.3.4.

The density adaptor, α, is defined as the density adaptor for the landslide mass and is the ratio of the accumulated mass density to the removed mass density. The mass density is related to the structure and composition of the soil and rock in the landslide area, which vary according to the location and the type of landslide. The value of the density adaptor is, therefore, generally obtained from field data.

In MBM, the calculation terminates when $V_k^{-}$ is equals to αV⁺, a theoretical condition achieved as long as the gap is smaller than a predefined threshold. The definition of an appropriate threshold is, therefore, crucial to the degree of convergence and accuracy. In addition, as explained earlier, the step size for H^b determines both V⁻ and V⁺. This means that an appropriate threshold is difficult to define and an alternative solution is, therefore, introduced below.

2.2.3. Program Termination Condition

It can be inferred from MBM, and the above analysis, that each change in H^b produces a pair of values of V⁻ and V⁺, such that an increase in one of the pair produces a decrease in the other. Therefore, as a result of the dynamic movement of H^b, the plotted curves of V⁻ and V⁺ should intersect exactly at the balance point; that is where $V_k^{-}$ is equals to αV⁺. This intersection can, therefore, be used as a termination condition. As can be seen from Figure 2b, curve A₀A₁A₂ is V⁻ and B₀B₁B₂ is αV⁺. The coordinates on these curves can then be expressed as A₀(x₀, y₀¹), A₁(x₁, y₁¹), A₂(x₂, y₂¹), B₀(x₀, y₀²), B₁(x₁, y₁²), B₂(x₂, y₂²), and O(x, y) (the intersection point). When the intersection of the two curves occurs in the first time, the termination condition is achieved as the following:

$$\left(y_1^1-y_1^2\right)\left(y_2^1-y_2^2\right)<0$$

(7)

2.2.4. Determination of the Final Results (H^b, V⁻, and V⁺)

As discussed above, the intersection of the plots of V⁻ and αV⁺ represents a balance position where V⁻ equals αV⁺. Four known points can define two straight lines and the coordinates of the point of intersection of these lines can then be calculated. The four points adjacent to O are defined as A₁(x₁, y₁¹), A₂(x₂, y₂¹), B₁(x₁, y₁²), and B₂(x₂, y₂²). A₁A₂ defines Line 1 and B₁B₂ defines Line 2; O(x, y), the point of intersection of Line 1 and Line 2, can then be expressed as follows:

$$y=\left(y_2^2{y}_1^1-y_2^1{y}_1^2\right)/\left(y_2^2-y_1^2+y_1^1-y_2^1\right)$$

(8)

$$x=x_1+\left(x_2-x_1\right)\left(y-y_1^2\right)/\left(y_2^2-y_1^2\right)$$

(9)

If in one step of H^b (from A₁ to A₂), V⁻ and V⁺ change in a linear fashion, then x in O(x, y) is the final value of H^b and y is the calculated V⁻. The accumulated volume, V⁺, is then compared by y and α.

2.3. Calculation Workflow

The overall workflow of the proposed approach is illustrated by Figure 3, in which the input data, include remote sensing images, pre- and post-landslide DTM data, and the density adaptor. The only control parameter that needs to be set is the size of the H^b step. The output data are the removed and accumulated volumes of the landslide together with H^b. There are seven key steps in the approach.

(1)

Definition of the landslide-affected area. The area affected by the landslide is defined manually using remote sensing images.

(2)

Preprocessing of the DTM data. Pre- and post-landslide DTM data are examined; if there is a discrepancy in the sampling intervals, a resample is conducted with the bilinear method.

(3)

Co-registration of the two DTMs. The pre- and post-landslide DTM data are match-optimized to eliminate co-registration errors. The classic normalized cross-correlation (NCC) algorithm for image matching is applied [48–50]. The translation that maximizes the correlation between the two DTM was searched. Due to the effect of topographical changes on the result of the NCC, the area defined by (1) is used to mask the DTM. In this way, within the area affected by the landslide, where the elevations according to the pre- and post-landslide DTMs are different will be eliminated by the NCC matching. This approach helps to improve the matching accuracy of the NCC.

(4)

Definition of the DTM data within the landslide-affected area. After (2) and (3), the pre- and post-landslide DTMs should have the same sampling interval and accurate co-registration. Optimized DTM data are then collected from the landslide-affected area defined in (1).

(5)

Calculation of the volumes V⁻ and V⁺. V⁻ and V⁺ are obtained using MBM and are then compared to decide whether H^b is positive or negative. If V⁻ is larger than V⁺, H^b is positive; if not, then H^b is negative. The step size of H^b and the density adaptor must be set.

(6)

Define whether there is an intersection between the plotted values of V⁻ and V⁺. The calculation terminates when the first intersection occurs.

(7)

The four values of V⁻ and αV⁺ from the last two calculation steps are used to derive intersection coordinates according to Equations (8) and (9). The horizontal coordinate, x, is the height adaptor and the vertical coordinate y is V⁻, then V⁺ is calculated by y and α.

3. Case Study

3.1. Study Area and Data Collection

The study area is located in Xuankou, Wenchuan County, on the west side of the Mingjiang River (N30°59′49.7″–31°00′4.2″, E103°27′40.3″–103°27′59.9″). This typical mountain terrain has an average altitude of 1,000 m (minimum 750 m and maximum 1,550 m), with an average gradient of 24° (maximum 80°). The average temperature is 10–14.4 °C and the average annual precipitation is 1,400 mm (most of which falls in July and August). This area was devastated by landslides and debris flows in the 12 May 2008, Wenchuan Earthquake, which had a magnitude of Ms 8.0 and induced more than 56,000 landslides in steep mountainous terrain covering an area of about 41,750 km² [51,52].

To test the accuracy and accessibility of MBM, we chose the landslide located at 30°59′57.92″N 103°27′53.07″E, which is located behind Younian Village, alongside Guxigou, a branch of the Minjiang River, 5.3 km to the south of the earthquake epicenter. This is a typical landslide triggered by the tremendous energy of the active Yingchuan–Beichuan fault in the Longmen Mountain fault zone (Figures 4 and 5). The slide mass is basically composed of loose gravel with sand and stones; its back and side walls are in the shape of round-backed chair and its ligulate front is covered with vegetation, stones, and fallen surface soil. The removed zone and the accumulated zone are distinctive and typical (Figures 4 and 5).

For this area, high-resolution aerial images, pre- and post-landslide DTM data and field investigation data were acquired.

Aerial stereo images were acquired by the ADS40 aerial camera system on 16 May 2008, four days after the landslide. The aerial stereo images were processed by Pixel Factory [53] into a Digital Orthophoto Map (DOM) (Figure 5) and DTM (Figure 6b). The DOM and DTM were both in geotiff format with WGS-84 coordinates. For the same area, the ASTER GDEM (Figure 6a) made before the landslide occurred was acquired from NASA through the NASA Land Processes Distributed Active Archive Center (LP DAAC). The estimated accuracy of this GDEM was 20 m at a 95% confidence for vertical data and 30 m at a 95% confidence for horizontal data. The ASTER GDEM was in geotiff format with geographic latitude/longitude coordinates and a 1 arc-second (30 m) grid [54]. A comparison between the pre- and post-landslide DTMs (Figure 6) clearly shows the topographical changes caused by the landslide.

Table 1. Basic information about the Digital Orthophoto Map (DOM) and digital terrain models (DTMs) and their accuracy.

Data Type	Grid Spacing	Accuracy	Data Acquisition
DOM-post	0.6 m	±2 m	ADS40 system
DTM-post	3 m	Horizontal ±3 m; vertical ±5 m	From ADS40 stereo images by Photogrammetry
DTM-pre	∼30 m	Horizontal ±30 m; vertical ±20 m	ASTER GDEM

contains basic information about the DOM and the pre- and post-landslide DTMs and their accuracy.

Table 2. Basic statistics for the pre- and post-landslide DTMs (units: m).

DTM Data	Min	Max	Mean	Stdev	Max Elevation Difference
DTM-pre	797	1,411	970	130	614
DTM-post	773	1,384	930	130	611

contains basic statistics derived from the pre- and post-landslide DTMs and shows that there was a systematic error of around 30 m in the DTM elevation.

In September 2009, field investigations were carried out to acquire samples for density adaptor analysis. The sampling site is located in the alpine gorge area with high mountains and deep valleys. Meanwhile, the selected landslide is a rock one with a relative homogeneity of rock and soil material. Thus, we had difficulty in conducting field sampling and selected eight typical samples, which can be divided into two sets representing the pre- and post-landslide material composition respectively (see the locations shown in Figure 5). Samples i and ii were taken from the upper stable area of the slide and represent the pre-landslide surface soil density, whereas samples iii and iv were taken from the bedrock of the back wall and represent the pre-landslide gravel density. Samples I and II were taken from the collapsed slide mass and represent the post-landslide rock mass density, whereas samples III and IV were taken from the flow slide mass and represent the post-landslide surface density. The mass and volume for all the samples were measured and the density calculated,

Table 3. Mass, volume, and density of the field samples.

Pre-Landslide Samples				Post-Landslide Samples
Sample ID	Mass (g)	Volume (cm3)	Density (kg/m3)	Sample ID	Mass (g)	Volume (cm3)	Density (kg/m3)
i	11,273.0	4,029.0	2.8 × 10−3	I	64,017.0	23,658.4	2.7 × 10−3
ii	18,446.0	6,475.3	2.8 × 10−3	II	16,998.0	6,510.4	2.6 × 10−3
iii	69,896.0	23,197.0	3.0 × 10−3	III	74,762.0	27,769.5	2.7 × 10−3
iv	62,637.0	20,175.3	3.1 × 10−3	IV	88,167.0	32,776.8	2.7 × 10−3

lists the results. The average density of samples i–iv, 2.9 × 10⁻³ kg/m³, was taken as the pre-landslide material density and the average of samples I–IV, 2.7 × 10⁻³ kg/m³, was taken as the post-landslide density. According to Equation (4), the value of the density adaptor was 0.9094, which was determined by comparing the average density of the samples I–IV and that of i–iv.

3.2. Calculation and Results

The volumes for the Xuankou landslides were then calculated according to the procedure described in Section 2.3.

First, the extent of the landslide was carefully drawn manually using ADS40 ortho-images and saved as vector data (*.shp) (The yellow circles in Figure 5 show the main landslide-affected area.) Next, pre- and post-landslide DTM data were processed by re-sampling and then by co-registration. Before registration, the ASTER GDEM was converted to WGS84 coordinates, the same as used in the post-landslide DTM. Resampling from high-resolution to low-resolution images causes a loss of information in the DTM, whereas with resampling from low-resolution to high-resolution data there is no information loss. Therefore, the horizontal sampling interval for the post-landslide DTM remained unchanged while the pre-landslide DTM sampling interval was resampled with bilinear method from 30 m to 3 m.

The NCC algorithm was adopted for the plane registration. The high-accuracy post-landslide DTM was used as a reference. The pre-landslide DTM was moved back and forth to register it to the post-landslide DTM within a range of 800 m in the up-down and left-right directions, which was derived from a priori knowledge. To get a better correlation coefficient for the overlapping area, the landslide extent was used as a mask in order to subtract the landslide area. Each movement produced one value for the correlation between the pre- and post-landslide DTMs. The maximum coefficient was 0.991, which occurred at the position (45 m, 200 m), which means that the pre- and post-landslide DTMs had vertical and horizontal gaps of 45 m and 200 m, respectively.

After the resampling and co-registration, the pre- and post-landslide DTM data that had been optimized in this way, and the *.shp files representing the approximate landslide extent, were used to extract the effective coverage of the pre- and post-landslide DTM data. MBM was then used to calculate the removed and accumulated volumes. The following two parameters were then set.

Choice of H^b step. The definition of such a step is actually a compromise between calculation accuracy and efficiency. The smaller the step, the higher the accuracy but with inevitable loss of efficiency. Conversely, less demand for accuracy will help to increase the efficiency. In the present case, H^b step was set at 0.1 m.

Choice of the density adaptor. In this case, the density adaptor was set at 0.9094, which was the value derived from the field samples.

When all the data were ready and the control parameters had been set, it was possible to calculate V⁺ and V⁻ and to use the program to search for the balance line automatically so that the accumulated volume could be calculated. Figure 7 shows all the results calculated during this process. As H^b increases, the post-landslide DTM elevation also increases with the removed volume dropping and the accumulated volume rising. When H^b reaches 52.4 m, V⁻ and αV⁺ intersect for the first time and the calculation terminates. The four values of V⁻ and αV⁺ of the last two calculation steps (See

Table 4. The final two pairs of V⁻ and αV⁺ Values at Mass Balance Model (MBM) Termination.

Hb (m)	V− (m3)	αV+ (m3)
52.3	1.2209 × 106	1.2107 × 106
52.4	1.2146 × 106	1.2240 × 106

) are imported into Equations (8) and (9) and the final results are then calculated. The final value of H^b is 52.35 m. The removed volume is 1.2176 × 10⁶ m³, whereas the accumulated volume equals 1.3389 × 10⁶ m³.

In addition to MBM, HDM (Equation (1)) and AHDM [27] were also used to calculate the Xuankou landslide volumes. First, the accumulated volume was calculated using HDM: in this case, the post-landslide DTM elevation is always lower than that of the pre-landslide DTM due to the systematic error. For this landslide, V⁺ was zero and V⁻ was equal to 1.0807 × 10⁶ m³. The second method used an improved version of Equation (1), AHDM, in which five sub-areas (boxes numbered 0 to 4 in Figure 5) where there were no changes in elevation were chosen as the control zone for the DTM height adjustment. The difference between pre- and post-landslide DTMs for each of these five sub-areas was calculated and the average called meanID0-meanID4. The overall average for these five sub-areas was called meanID-ALL. All of these mean values were used as the height adjustment used to correct for the systematic difference between the two DTMs. V⁻ and V⁺ were then inferred when the systematic difference between the pre- and post-landslide DTMs was removed. All the results of HDM, AHDM, and MBM are summarized in

Table 5. Results obtained using different models and for different conditions.

Models		Height Adjustment (m)	V+ (×106 m3)	V − (×106 m3)
HDM		0	0	1.0807
AHDM	meanID0	52.5265	1.2068	1.3645
	meanID1	50.9075	1.3133	1.1330
	meanID2	53.0371	1.1759	1.4402
	meanID3	40.6712	2.4804	1.6327
	meanID4	30.3309	4.4775	1.7860
	meanIDAll	45.4946	1.7977	4.8748
MBM		52.3524	1.2176	1.3389

.

3.3. Discussion

3.3.1. Co-Registration of Pre- and Post-Landslide DTMs

In both the conventional (including HDM and AHDM) and the improved model (MBM), the elevation difference is an essential element for the accumulated volume calculation and the results are extremely sensitive to errors in the co-registration. In fact, pre- and post-landslide DTMs may differ in terms of the coordinate system, positional accuracy, the baseline of the elevation, and other aspects. Therefore, co-registration is indispensable. In order to demonstrate this issue, the elevation differences between the pre-and post-landslide DTMs before and after the co-registration were calculated and the histograms of the elevation difference images were also calculated. The results are shown in Figure 8. The histogram before registration has two peaks and the overall values range from −234.95 to 186.23; the histogram after registration has only one peak and the values range from −162.31 to 48.82 with most of the values being found in the range (−75, 0), which is clearly more concentrated. Clearly, the distribution of the latter histogram is much more reasonable than that of the former. These observations verify the significance and effect of co-registration. In addition, in Figure 8, the peak in the post-registration histogram appears around −26, not the theoretical zero position, which means that that there was a systematic elevation difference between the pre- and post-landslide DTMs of around 26 m.

3.3.2. Reliability and Accuracy: Improved Model against Conventional Model

As shown in Table 5, due to height value (Z-coordinate) errors in the pre- and post-landslide DTMs, without optimization, HDM cannot give the accumulated volume and the calculated removed volume has no significance. The results show that only when the input data are ideal does it make sense to apply this equation and, thus, its application is very limited.

When AHDM is applied, the value of the height adjustment is obtained from the average values for selected sub-areas. As can be seen from Table 5, when a different sub-area is selected, a different height adjustment value can be obtained and V⁺ and V⁻ will also change significantly. This means that AHDM has poor reliability and low accuracy. Only when the correct height adjustment value is used can AHDM produce good results. In fact, the elevation differences between the pre- and post-landslide DTMs may range from small to large or even huge, indicating the impossibility of choosing ideal elevation control sub-areas.

In contrast to AHDM, the height adjustment in MBM is automatically derived with the constraints that the amount of removed material equals the amount of accumulated material. It has a unique value and, as long as the pre- and post-landslide DTMs have no errors in plane position, H^b should be correct. In this sense, MBM eliminates the error due to the elevation difference between the pre- and post-landslide DTMs and shows good reliability. In addition, V⁺ and V⁻ are calculated by the same method using the same data, which means that the sources of possible error for V⁺ and V⁻ are almost the same. In other words, the errors in V⁺ and V⁻, which are present on both sides of Equation (5), are of approximate size. This means that it is now possible to find the mass balance line correctly when the input data contain errors. It can, therefore, also be inferred that the constraint of the mass balance improves the reliability and accuracy of MBM in this way.

Figure 9 is the profile corresponding to the pink dotted line in Figure 5. It demonstrates the advantage of using MBM instead of the conventional model. Initially, the post-landslide DTM hatching line follows the blue curve. Due to the elevation difference errors, the post-landslide DTM profile always lies below that of the pre-landslide DTM. When HDM is applied here, no accumulated volume can be obtained and so V⁺ is 0. When AHDM is applied, 45.49 m is added to the post-landslide DTM profile (and meanID-ALL is used as the height adjustment value). Now represented by the green line, this new profile intersects the pre-landslide curve at O₂. When MBM is applied, 52.35 m is added to the post-landslide DTM profile and the new profile, represented by the red line, intersects the pre-landslide curve at O₁.

3.3.3. Evidence for the Mass Balance Line

The mass balance line was defined as the boundary at which the mass of the landslide-triggered accumulated material equals that of the material removed from the upper slope. Another possible definition is that this is the line where the elevations in the pre- and post-landslide DTMs are equal. This means the mass balance line lies along the contour line with zero elevation difference values. The mass balance line is a theoretical idea and cannot be found by field-work; however, it can be found from Figure 5.

The red dotted curve in Figure 5 is the contour line for zero elevation difference according to AHDM (with mean ID-ALL used as the height adjustment value). Within this contour line the elevation difference is larger than zero, representing the accumulated area. The area outside the line is the removed area, where the elevation difference is smaller than zero. The red curve in Figure 5 is the contour line marking zero elevation difference according to MBM, inside this line the elevation difference is smaller than zero, representing the removed area. The area outside represents the accumulated area, where the elevation difference is larger than 0. According to the definition of the mass balance line, the lower boundary of zero elevation difference contour according to MBM is the mass balance line of the study case, which is blue thick line in Figure 5. In other words, the reality of the situation is entirely consistent with the theory of the proposed model. This evidence for the mass balance line supports the concept of the mass balance principle and the reliability of MBM.

A comparison of the two contour lines in Figure 5 very clearly shows that the result of applying MBM is much more reasonable than that obtained using AHDM. When the former is applied, the removed area is in the upper part of the landslide and lies along the landslide boundary. In addition, the area of removed material is smaller than that of accumulated material, which is consistent with Figure 9. We can also see from Figure 9 that above the intersection point (in the range between A and O₁), the elevation differences tend to be larger, while, below the intersection point, they tend to be smaller. One possible reason is that when the landslide occurs, the mountain collapses and the lower part of the deposit is more dispersed and loose than the upper part. The 3D image (Figure 4) supports this idea.

3.3.4. Parameter Setting in MBM

(a) Setting of H^b

H^b, the driving force behind MBM’s automatic search for the mass balance line, is a key parameter in that equation. Figure 7 shows that V⁻ and V⁺ change simultaneously and finally achieve balance with H^b. The dynamic change in the values of V⁻ and V⁺ stands as solid proof of the objectivity of the mass balance principle. The automatic search for the balance line is another advantage of MBM. Once the model is set up, the program can automatically define a direction of movement with V⁻ and V⁺ being calculated accordingly.

In MBM, the only parameter to be determined is the step size for H^b. Clearly, an increase in this interval means higher efficiency and lower accuracy. However, if the interval is too large then the assumption contained in Equations (8) and (9)—that V⁻ and V⁺ show linear changes for one step—is invalid. To investigate the appropriate size of the H^b step, the results of calculations with the H^b step ranging from 5 mm to 1 m were compared (

Table 6. The computational accuracy and efficiency of MBM with different H^b step value.

Step (m)	0.005	0.01	0.05	0.1	0.3	0.5	1
Hb (m)	52.3524	52.3524	52.3524	52.3524	52.3524	52.3524	52.3526
V+ (×106 m3)	1.3389	1.3389	1.3389	1.3389	1.3390	1.3391	1.3395
Number of loops	10,471	5,236	1,048	524	175	105	53
Computation time (second)	2,241	1,120	227	103	39	24	13

Note: The computation was processed in Matlab (Version R2011a) by personal PC (Processor: Intel^® Core i5-2400 [email protected] GHz. RAM: 4.00 GB).

). As the H^b step increased 200 times from 5 mm to 1 m, the final H^b value increased by 0.2 mm and V⁺ increased by 600 m³, representing changes of 0.0004% and 0.04%, respectively. These subtle changes indicate that the final results of MBM (H^b, V⁻, and V⁺) are not sensitive to the size of the H^b step. As to the computational efficiency, the smaller the H^b step, the bigger computational cost: Table 6 gives the computation time and the number of search loops used to find the balance line.

It can, therefore, be inferred that, under the balance condition, the assumptions of Equations (8) and (9) are valid. In practice, it is recommended that the size of the H^b step be set between 0.005 and 0.1 m so that computational efficiency and accuracy are given equal importance.

(b) Setting of density adaptor

The density adaptor is the ratio of the accumulated mass density to the removed mass density and is used to adjust for the difference in the mass density before and after the landslide. As mentioned earlier, the value of the density adaptor varies according to the material type and also the composition and structure of the landslide rock and soil. This value is generally obtained from field data. In this study, we mainly gave an illustration of the methods used to take samples for the calculation of the density adaptor and also the methods used to calculate it.

From the role of the density adaptor in Equation (5) and the pseudo-code 6, the density adaptor can also be served as a regulation parameter of the volume estimations in MBM. When MBM is used to estimate the volume of a landslide where the deposit has been carried away and the material balance destroyed, the density adaptor can also be used to adjust for the removed mass. If the ratio of the lost mass to the removed mass is estimated through field investigation and survey, a new mass balance can be derived as Equation (10). In this way, the error due to the lost mass can be reduced to a certain degree and MBM can be made suitable for application to a wider range of landslides.

$$V^{-}=\frac{a}{(1-b)}{V}^{+}=a^{\prime }{V}^{+}$$

(10)

where, a is the density adaptor, b is the ratio of the lost mass to the accumulated mass, and α′ can be seen as the new density adaptor, which takes the loss mass into account.

3.3.5. Limitations and Errors of MBM

As with the conventional model, MBM cannot be applied to all situations. Strictly speaking, it applies only to landslides with no material loss and is more applicable to rotational and translational slides of more solid materials. Where the deposit is carried away and the material balance is destroyed, MBM will produce nothing but a smaller accumulated volume. By field investigation and survey, a new density adaptor may be estimated, and the error due to the lost mass can be reduced to a certain degree.

In addition, the result produced by MBM is not the overall volume of the deposited material (total deposited volume in Figure 1a) but only the volume of the region where surface morphology deformation occurs and the elevation before and after the landslide event is different. As for the “filling up” volume, although it is part of the deposit, this volume cannot be calculated using MBM because there are no differences in the elevation before and after the landslide. MBM cannot, therefore, be used to calculate the total deposited volume for the landslide.

Apart from the limitations of MBM, there are some errors related to the input data, the remote sensing imagery used and the pre- and post-landslide DTMs. The mechanism by which the uncertainties in the remote sensing and DTM data propagate through to the error in the final volume calculation is very complex and only a brief analysis of the error propagation is provided here.

The role of the image is to define the landslide-affected area, which then is used to select the pre- and post-landslide DTMs. Many factors related to the image, such as the resolution, the plane position, and the co-registration with the DTMs, will affect the range of the selected DTMs and lead to a reduction in the accuracy of the final results. To reduce the error related to the definition of the landslide-affected area, a larger area than the true affected area can be defined. As the difference between the pre- and post-landslide DTMs outside the affected area is almost zero, MBM is less sensitive to the input data used, as mentioned in Section 3.3.2.

The uncertainty in the DTM data is due to many factors, such as the original data used, the model, the interpolation algorithm, the mesh grid distance, the height baseline. This uncertainty affects the accuracy of the final results (the removed volume and the accumulated volume) in two respects. Firstly, the uncertainty reduces the authenticity of the morphology that was built by the DTM data. The gap between the real morphology and the DTM data cannot be avoided or eliminated by either the conventional model (including HDM and AHDM) or MBM. Secondly, the uncertainty increases the inconsistency between the different versions of the pre- and post-landslide DTMs and then reduces the accuracy of the final landslide volume results. In order to control the uncertainty in the final results, the user of MBM is recommended to select DTMs that have a high degree of consistency regarding the plane position, mesh grid distance, height baseline, etc.

Differences in the DTM sampling interval is a common inconsistency that should be taken into account when different versions of DTMs are used to calculate the landslide volume. Resampling from high-resolution to low-resolution data causes a loss of information in the DTM, whereas, with resampling from low-resolution to high-resolution data, there is no information loss. The latter method is, therefore, recommended when processing different versions of DTMs.

4. Conclusions

In this study, the mechanisms and laws governing the material deposited by landslides are studied and then the mass balance principle and the mass balance line are defined. By the mass balance principle, the mass of the accumulated material is the same as that of the removed material on the upper slope before the landslide event; the mass balance line is defined as the boundary at which the mass of the accumulated material equals that of the removed material. Based on the mass balance principle and the definition of the mass balance line, a novel and improved model (MBM) for calculating the accumulated volumes of landslides is proposed. Due to the use of mass balance constraints, the height adaptor, and the volume difference as grounds for program termination, MBM can eliminate or reduce some of the errors that occur in the conventional methodology, such as the differences in data accuracy and the difference in the height baseline of the pre- and post-landslide DTMs. MBM is, therefore, much more reliable and is simple to implement. By applying a parameter named the height adaptor, MBM translates the volume calculation into an automatic search for the mass balance line within the scope of the landslide.

However, as for the conventional model, MBM cannot be applied in all cases and applies well only to landslides where there is no material loss. If, instead, some of the deposit is carried away and the material balance is destroyed, more field investigation is needed to carry out to produce the new density adaptor according to Equation (10). In addition, MBM can only calculate the volume of the region where the surface morphology deformation occurs and the elevation is different before and after the landslide event (the accumulated volume in Figure 1a). As for the “filling up” volume (see Figure 1a), although it is part of the deposit, as there are no differences in the elevation before and after the landslide, this volume cannot be calculated.