AscDAMs: advanced SLAM-based channel detection and mapping system

Obtaining high-resolution, accurate channel topography and deposit conditions has been a challenge for the study of channelized debris flow. Currently, widely used mapping technologies including satellite imaging and drone photogrammetry struggle to precisely observe channel interior conditions of long and deep mountainous gullies, particularly those in the Wenchuan earthquake region. SLAM is an emerging tech for 3D mapping; however, extremely rugged environment in long and deep gullies poses two major challenges even for the state-of-the-art SLAM: (1) atypical features and (2) violent swaying and oscillation of sensors. These issues result in large deviation and lots of noise for SLAM results. To improve SLAM mapping in such environments, we propose an advanced SLAM-based channel detection and mapping system, namely AscDAMs. It features three main enhancements to post-process SLAM results: (1) the digital orthophoto map-aided deviation correction algorithm greatly eliminates the systematic error; (2) the point cloud smoothing algorithm substantially diminishes noise; (3) the cross-section extraction algorithm enables the quantitative assessment of channel deposits and their changes. Two field experiments were conducted in Chutou gully, Wenchuan County in China in February and November 2023, representing observations before and after the rainy season. We demonstrate the capability of AscDAMs to greatly improve SLAM results, promoting SLAM for mapping the specially challenging environment. The proposed method compensates for the insufficiencies of existing technologies in detecting debris flow channel interiors including detailed channel morphology, erosion patterns, deposit distinction, volume estimation and change detection. It serves to enhance the study of full-scale debris flow mechanisms, long-term post-seismic evolution, and hazard assessment.

1 Introduction

The 2008 Wenchuan earthquake produced a huge amount of loose solid material, spawning repeated post-seismic channelized debris flows in numerous long and deep gullies (Tang et al., 2009; Guo et al., 2016; Fan et al., 2019), which have posed a continuing threat to human lives and properties (Xu et al., 2012; Hu and Huang, 2017; Zhang and Zhang, 2017; Fan et al., 2018, 2019; Shen et al., 2020; X. Z. Zhang et al., 2022; Zhang et al., 2023). During the 15 years after the earthquake, the loose hillslope material had been transported gradually into channel deposit (Zhang and Zhang, 2017), resulting in significant change in the occurrence frequency and initiation mechanisms of debris flow (Berger et al., 2011b; Chen et al., 2024). The channel morphology continually changes with the process of sediment transportation. The morphology of debris flow channel indicates the channel topography elevation, loose material distribution, debris flow impact area, entrainment area and entrainment depth, deposition area and volume, etc. (Remaître et al., 2005). An accurate observation of debris flow channel morphology is vital to the deep understanding of the debris flow initiation mechanism and risk assessment. Besides, the accuracy of many debris flow numerical simulation frameworks, e.g. r.avaflow (Mergili et al., 2017) and EDDA (Chen and Zhang, 2015; Shen et al., 2018), and debris flow assessment factors (Liang et al., 2012; Meyer et al., 2014; Li et al., 2021) relies on the accuracy of topography data.

Conventionally, morphology data of debris flow channels are mainly obtained through satellite images and field investigation. Satellites can easily produce digital elevation models (DEMs) and digital orthophoto maps (DOMs) (Zhang et al., 2014; Mueting et al., 2021; Luo et al., 2022) of a wide range of area (over tens of square kilometres) with low average unit cost. However, satellite-derived DEMs and DOMs often exhibit limited resolution and accuracy in rugged terrains, especially in mountainous areas and deep valleys (Sun et al., 2015; Zhou et al., 2017; Liu et al., 2021). Field investigation by visual observation and manual measurement by ruler or laser range finder could be less efficient and imprecise. Recently, equipment like radars (Schurch et al., 2011; Caduff et al., 2015; Morino et al., 2019; Bonneau et al., 2022) and UAVs (Simoni et al., 2020; Walter et al., 2022; X. Z. Zhang et al., 2022; Zhang et al., 2023) has been widely applied in field investigations to get more detailed topographic maps. However, these methods are still restricted by the application environment. For example, UAV mapping requires open airspace, enough signal strength of GNSS, and skilled operators (Cucchiaro et al., 2019; Imaizumi et al., 2019; Huang et al., 2022). It is extremely difficult to operate in narrowly spaced, signal-blocked, and GNSS-denied deep valley environments. The radar-based approaches require an appropriate arrangement for installation locations and/or scanned areas, which is also very difficult in this rugged and rocky region (Blasone et al., 2014; Caduff et al., 2015; Morino et al., 2019; Tang et al., 2022). Hence, current technologies could not provide sufficient, accurate, and consistent information in deep valleys of alpine areas, where channelized debris flows initiated and developed, especially the Wenchuan earthquake region. These constraints impede a thorough understanding of the debris flow mechanisms. Therefore, developing a new method for accurately detecting debris flow channels is an urgent, common key issue in channelized debris flow research and hazard mitigation.

Simultaneous localization and mapping (SLAM) (Bailey and Durrant-Whyte, 2006; Durrant-Whyte and Bailey, 2006; Cadena et al., 2016; Barros et al., 2022) technology is a mobile measurement method that continuously records data on the move. It has wide usage in the field of robotic navigation, autonomous driving, and topographic mapping. One typical kind of SLAM technology is lidar odometry and mapping (LOAM) (Zhang and Singh, 2014) based on light radar (lidar), which matches point clouds with their features. LOAM distinguishes features based on the curvature of points scanned by lidar. By this method, the computational complexity can be reduced. SLAM technology has been applied across different platforms for numerous scenarios of topographic mapping, including forestry (Kukko et al., 2017; Pierzchala et al., 2018; Li et al., 2023), underground tunnels (Ullman et al., 2023), urban morphology (Tanduo et al., 2022), hillslope gullies (Kinsey-Henderson et al., 2021), and densely vegetated hillsides (Marotta et al., 2021). SLAM has excellent potential to supplement the data of debris flow channel conditions where the satellite images are of low quality.

Nonetheless, it is still challenging to apply state-of-the-art SLAM techniques in debris flow gullies in alpine areas (e.g. Wenchuan earthquake region) from in situ preliminary tests. Long and deep debris flow gullies with extremely rugged environments have two major challenges for successful application of SLAM (see visual details in Sect. 2). One is atypical features. The channelized debris flow gullies are typically long and deep. Channel sidewalls are sparsely vegetated, while the channel bed is bestrewed with loose materials. The features of this specific environment are atypical. It is difficult to acquire sufficient and effective channel morphology data. The second is violent swaying and oscillation of sensors. Due to the presence of large rocks and flowing streams within the channel, activities like climbing, jumping, and rotating are unavoidable during the data collection process. Hence, the sensors experience inevitable severe swaying and rocking. These issues result in the following: firstly, the SLAM algorithm was apt to producing large mapping deviation because the algorithm is unable to extract enough effective features for matching and computing in this environment with atypical features. Besides, SLAM algorithms relying on scanning and matching of point cloud frames will lead to systematic tiny errors that accumulate into considerable deviation with the increase of channel length and the expanding of the dataset. Although many efforts have been made to mitigate the influence of the deviation, for example, introducing visual–inertial odometry, looper closure, and GNSS information for pose correction (Shan et al., 2020; Lin and Zhang, 2022), the mapping deviation is still prominent for such long and deep gullies. Secondly, due to the significant oscillation of sensors on the backpack-type collection system during mountaineering process, the mapping result of SLAM algorithm contains lots of noise.

The above limitations have hindered the application of current SLAM technology in investigating debris flow channels. To address these problems, we propose an advanced SLAM-based channel detection and mapping system (AscDAMs), which contains three major novel contributions to post-process SLAM results.

Figure 1 A typical mountainous area in Wenchuan County, Sichuan, China: (a) 3D map of the mountainous area of Wenchuan on 14 November 2021 from © Google Earth, (b) overview of Chutou gully on 25 March 2020 from ZY03 satellite image, and (c) accumulation zone on 14 November 2021 from © Google Earth.

1. Deviation correction algorithm. The newly proposed algorithm uses the DOM and barometric elevation data as a benchmark to minimize the accumulation drift, offering a viable solution to the morphology detection challenge in long and deep channels.

2. Point cloud smoothing algorithm. This algorithm is designed to mitigate the noise caused by drastic motion state change. It has the capability to enhance the quality of SLAM results, even in scenarios where sensors undergo violent oscillations during data detection in the harsh channel environment.

3. Cross-section extraction algorithm. This is beneficial for extracting typical channel cross-section profiles. This functionality enables the quantitative assessment of the channel interior, including the channel deposit distinction, volume estimation, change monitoring, and erosion pattern observation.

In the following content, the algorithms and application process of AscDAMs will be introduced in detail. The outputs of AscDAMs are demonstrated to have huge potential to greatly propel the comprehensive exploration of the underlying mechanisms of debris flow and the long-term evolution of debris flow activities.

Figure 2 Environment inside the Chutou gully on 10 February 2023.

2 Application environment

As one of the typical post-seismic debris flow gullies in Wenchuan County, the Chutou gully is featured by steep hillslopes and a deep channel (see Fig. 1a and b). The disastrous earthquake on 12 May 2008 provided a huge amount of loose material inside the gully that can be easily mobilized into debris flows (Tang et al., 2009; Guo et al., 2016; Yang et al., 2023a). Three large-scale debris flows occurred in the Chutou gully catchment during the rainy season in 2013, 2019, and 2020 (X. Z. Zhang et al., 2022; Zhang et al., 2023). The repeated debris flow hazards not only pose a continuing threat to local people's lives and properties but also change the channel's and Min River's morphology completely (see Fig. 1c). After the earthquake, most of the loose material on the hillslopes had moved down into the channel (Xiong et al., 2022; Y. Y. Zhang et al., 2022; Zhang et al., 2023). A recent conceptual model of the post-seismic debris flow shows that the debris flow activity in the Wenchuan earthquake region is still at the unstable stage, showing a relatively high future risk of debris flow (Yang et al., 2023b). This makes it imperative to know the channel interior conditions for mechanism study, debris flow prediction, and risk assessment.

The gully can be divided into a formation zone, a circulation zone, and a deposit fan (Fig. 1b). The formation zone is about 16 km² with complicated topographical conditions. We further partition the formation zone into a converging zone and an initiation zone. The converging zone is mainly exposed hillslopes for runoff convergence, which can be detected directly by satellites. Recent debris flows, occurring in 2019 and 2020, originated in the initiation zone (X. Z. Zhang et al., 2022; Zhang et al., 2023). Loose materials distributed in this zone play a pivotal role in initiating and amplifying debris flows. Downstream, the channel in the circulation zone is broad, and the terrain in the deposit fan is gentle, making detection easy by satellites or UAVs. However, different from these zones, the initiation zone is characterized by high-rising, inward-sloping, overhanging channel side walls, as indicated by Fig. 2. The overhanging cliffs obstruct the overhead view, and GNSS signal blockage by high mountains makes it difficult to obtain an accurate and detailed channel interior morphology by satellites and UAVs. Hence, we selected the initiation zone as the study area for the current research to showcase the novelty of the proposed approach.

SLAM algorithms were mainly applied to relatively stable environments (e.g. urban roads, river courses, park trails). They perform satisfactorily in these smooth environments (Ye et al., 2019; Shan et al., 2020). However, as depicted in Fig. 2, the channel bed is covered with loose deposits, rendering the channel path an exceedingly rugged terrain. The rocky cliffs on both sides of the channel are uneven. A natural stream exists inside the channel, but the water surface is also unsteady due to the channel gradient and bumpy bed. These atypical channel features make it particularly harsh for the implementation of SLAM algorithms. The long and deep channel morphology, together with the complex channel environment, extends the working duration of SLAM platforms (e.g. handheld, backpack, helmet), consequently enlarging the size of the dataset. This, in turn, will introduce additional errors to the results of SLAM. During the data acquisition, the huge stones and stream have emerged as obstacles, introducing difficulties to SLAM detection. The abrupt changes in the sensor pose caused by climbing and jumping, coupled with sudden swerving, can result in inaccurate pose estimations, consequently generating a significant amount of noise. AscDAMs is thus proposed to conquer these technical challenges for SLAM.

3 Methodology

3.1 Equipment and data

The data acquisition system (Fig. 3) of AscDAMs is a backpack type and consists of a lidar, an inertial measuring unit (IMU), a small-form-factor computer, a camera, etc. It is cost-effective and lightweight. The total weight of the device together with the battery is less than 10 kg, with the main cost not exceeding USD 10 000. Detailed information about its core configuration is shown in Table 1. We choose a lidar with 16 scan lines, an IMU with a frequency of 400 Hz, and a camera with a 1280×720 resolution. When assembling, it is important to note that the lidar and the IMU must be fixed on the bracket to avoid unexpected errors. The view of the lidar and camera should be unobstructed. The lidar and the IMU offer data for SLAM to calculate the map, and the camera is designed to provide colour information for channel characteristic extraction (e.g. channel deposit, vegetation). We utilized the system to record data of the selected zone on 10 February 2023, with the recorded section more than 1 km long (shown in Fig. 1b). A repeated channel scanning was conducted on 8 November 2023. It is noted that the second investigation is shorter than the first one because we were blocked by the rising water. The data characteristics are shown in Table 2.

Figure 3 The photo of the data acquisition system of AscDAMs.

Table 1 The list of core devices and their parameters.

Table 2 The characteristics of the recorded data in Chutou gully.

3.2 Algorithms

The entire workflow (Fig. 4) of AscDAMs consists of the application of SLAM and several newly proposed algorithms including deviation correction, point cloud smoothing, cross-section extraction, and map colouring. First, SLAM registers and fuses lidar point cloud frames with motion state (i.e. the acceleration and pose of lidar) recorded by the IMU to generate an initial global map and an initial trajectory of data acquisition. Then, we introduced a new algorithm for deviation correction of the initial global map, called the DOM- and barometer-aided deviation correction (DBADC). The coordinate of each point in the global map and trajectory will be corrected by referring to the horizontal coordinate of the DOM and elevation of barometer. Next, the modified global map will be smoothed utilizing a new algorithm, called the weighted-elastic-voxel grid (WEVG). Finally, cross-section extraction is processed, including cutting and projecting the smoothed global map and smoothing and densifying cross-sections. The final global map containing coordinate information can be converted to DEMs. In addition, a map colouring algorithm is developed to supplement more information and the global map. It fuses camera RGB data with lidar data and expands SLAM to obtain coloured maps. The final-coloured maps are obtained after colour optimization.

Figure 4 Algorithm flow scheme of AscDAMs. The initial SLAM algorithm calculates the initial maps utilizing the IMU and lidar. After optimizing the maps by deviation correction, point cloud smoothing, cross-section extraction, and map colouring, the channel global map with and without RGB colour can be obtained and employed in debris flow study.

3.2.1 Application of LIOSAM and map colouring

As one of the advanced LOAM algorithms, LIOSAM (Shan et al., 2020) merges with multiple sensors and enjoys the benefits of being accurate, open-source, and easy to replicate, which was widely used in 3D mapping. Therefore, this study is based on LIOSAM. Nonetheless, the choice of SLAM algorithm is not exclusive. AscDAMs is a universal channel morphology data acquisition and processing system, which is open to other LOAM algorithms. The success of applying LIOSAM to calculate the initial global map and trajectory map is a vital premise, as the subsequent optimization focusses on the two maps. For rugged debris flow channels, after various preliminary tests, it was found that the horizontal resolution of the input lidar point cloud images and the loop closure parameters have the greatest impact on whether we can obtain the maps. We achieved the application of LIOSAM in rugged environments by downsampling the lidar point cloud and employing loop closure as an iterative closest point (ICP) to rematch point clouds.

The initial LIOSAM map contains no colour information, while other SLAM algorithms have complex operations and strict requirements of equipment (Lin and Zhang, 2022; Zheng et al., 2022). Hence, we propose a procedure to get coloured maps of the channel. The procedure consists of three steps: colour fusion, RGB–SLAM, and colour optimization.

The method of colour fusion refers to an existing algorithm of fusing a camera and lidar on https://github.com/KAI-yq/camera_lidar_fusion (last access: 11 September 2024). We modify its way of colouring from point by point to frame by frame. This algorithm can fuse each lidar point cloud frame with a simultaneous camera picture to get a single-frame coloured point cloud, by calibrating and synchronizing the camera and lidar. RGB–SLAM is newly proposed to calculate the coloured map in the present study. It is developed by modifying the data acquisition interface of LIOSAM. The implementation of RGB–SLAM is similar to that of LIOSAM.

However, there are many uncoloured points in the initial map $C=\mathit{\{}{c}_{\mathrm{1}},c_{\mathrm{2}},c_{\mathrm{3}},\mathrm{\cdots },c_i,\mathrm{\cdots },c_I\mathit{\}}$ calculated by RGB–SLAM because of the limitation of the camera view. Hence, we have developed an algorithm to optimize the colour of the map. First, the uncoloured points in C are filtered to produce a new coloured map $C^{\prime }=\mathit{\{}{c}_{\mathrm{1}}^{\prime },c_{\mathrm{2}}^{\prime },c_{\mathrm{3}}^{\prime },\mathrm{\cdots },c_j^{\prime },\mathrm{\cdots },c_J^{\prime }\mathit{\}}$ as a colouring reference. Then, the colour data around the uncoloured point are searched out employing the K-nearest neighbour (KNN) method. The searched colour is given to the untouched point. Finally, points that still do not get colour after re-colouring are filtered out, and the final colourful map $C^{\prime \prime }=\mathit{\{}{c}_{\mathrm{1}}^{\prime \prime },c_{\mathrm{2}}^{\prime \prime },c_{\mathrm{3}}^{\prime \prime },\mathrm{\cdots },c_k^{\prime \prime },\mathrm{\cdots },c_K^{\prime \prime }\mathit{\}}$ is completed.

3.2.2 Deviation correction

The initial maps produced by LIOSAM have large deviations due to the accumulation of atypical features matching errors and systematical errors. These errors can cause the SLAM-generated map to drift from its actual position. Hence, it is necessary to correct the initial mapping to mitigate the drift. We introduced a new algorithm to correct the mapping deviation referring to the horizontal coordinate of the DOM and elevation of the barometer, called the DOM- and barometer-aided deviation correction (DBADC).

Theoretically, the initial maps, which contain stereo spatial data, need to be modified in X–Y–Z dimensions. However, for rotating lidars, the errors in the horizontal (i.e. X–Y) directions are the same and can be adjusted together. Hence, the proposed deviation correction algorithm introduces two scaling factors named horizontal scaling factor f_h and elevation scaling factor f_e to rectify the drift in the initial maps. Since GNSS data are not available in such gully environment, we estimate f_h with the satellite DOM and calculate f_e by converting the air pressure values recorded by the barometer embedded in the IMU to altitude. The scaling factors indicate the degree of modification between the initial point cloud map and the real channel morphology data. The closer the factors are to 1, the smaller the adjustment will be.

We develop an auto-estimation algorithm to determine f_h. Firstly, an adaptive segmentation method (Bradley and Roth, 2007) is employed to automatically extract debris flow channel area Ω from the DOM. Then the channel can be fitted as a line T_DOM by removing noise and curve fitting. After comparing every fragment of T_DOM with the initial trajectory $T=\mathit{\{}{t}_{\mathrm{1}},t_{\mathrm{2}},t_{\mathrm{3}},\mathrm{\cdots },t_i,\mathrm{\cdots },t_I\mathit{\}}$ and calculating their correlation, we can find out the most relevant fragment $T_{\mathrm{DOM}}^{\prime }$ to T (Eq. 1). $T_{\mathrm{DOM}}^{\prime }$ can be regarded as the correction benchmark of T. Finally, we utilize the ICP method to determine f_h (Eq. 2). After scaling the trajectory T with ratio variable s, this makes the difference between the scaled trajectory and $T_{\mathrm{DOM}}^{\prime }$ the smallest equal to f_h. The elevation scaling factor f_e can be determined by referring to the elevation of two real locations where the data acquisition system starts and stops recording (named “referred start point” and “referred end point” respectively), shown as Eq. (3). ΔL_e represents the elevation difference between these two points in the initial trajectory of the point cloud map, while L_e is their real elevation variation converted by barometer data.

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{T}_{\mathrm{DOM}}^{\prime }=\underset{T_{\mathrm{DOM}\mathrm{fragment}}\in T_{\mathrm{DOM}}}{\mathrm{argmax}}& \left[\left(T^{*}\cdot T_{\mathrm{DOM}\mathrm{fragment}}\right)/\right.\\ & \left.\left(∥T∥\cdot ∥T_{\mathrm{DOM}\mathrm{fragment}}∥\right)\right]\end{array}\text{(2)}& {\displaystyle }{f}_{\mathrm{h}}=\underset{s}{\mathrm{argmin}}\left[sT-T_{\mathrm{DOM}}^{\prime }\right]\text{(3)}& {\displaystyle }{f}_{\mathrm{e}}=L_{\mathrm{e}}/\mathrm{\Delta }{L}_{\mathrm{e}}\end{array}$$

In deviation correction of the initial global map, we group the points within the same spatial interval on the initial global map into the same correction unit. The coordinate deviations of all points within the same unit are considered to be the same. The shorter the interval length, the more accurate the correction result. There are two hypotheses for the units:

1. The construction and size of each unit are accurate. All the points in the same unit should be regarded as a whole and translated together when rectifying.

2. The overall errors of the initial global map are evenly distributed across the units.

The rationality of the two hypotheses is that we have maximized the performance of LIOSAM by downsampling and employing loop closure as enhanced ICP registration (see Sect. 3.2.1). The interval length of the unit is short enough so that its error could be neglected. However, the cumulative deviation by numerous units should be treated seriously.

We first smooth the trajectory T by employing the Fourier transform (and/or other smooth methods) $\widehat{f}$ to remove jitter. Then it is densified by linear interpolation $\widehat{L}$ to obtain the smoothed–densified trajectory $T^{\prime }=\mathit{\{}{t}_{\mathrm{1}}^{\prime },t_{\mathrm{2}}^{\prime },t_{\mathrm{3}}^{\prime },\mathrm{\cdots },t_j^{\prime },\mathrm{\cdots },t_J^{\prime }\mathit{\}}$ (Eq. 4).

$$\begin{array}{}\text{(4)}& T^{\prime }=\widehat{L}\left(\widehat{f}\left[T\right]\right)\end{array}$$

After projecting T^′ onto the horizontal (X–Y) plane, we get new trajectory point cloud $T^{\prime \prime }=\mathit{\{}{t}_{\mathrm{1}}^{\prime \prime },t_{\mathrm{2}}^{\prime \prime },t_{\mathrm{3}}^{\prime \prime },\mathrm{\cdots },t_j^{\prime \prime },\mathrm{\cdots },t_J^{\prime \prime }\mathit{\}}$ (Eq. 5), which can be considered as a division index trajectory.

$$\begin{array}{}\text{(5)}& t_j^{\prime \prime }=t_j^{\prime }\cdot \left[\begin{array}{ccc}\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right].\end{array}$$

The division of the initial global map is the key step of deviation correction. The adjacent points should be grouped into the same or adjacent unit. Next, we employ the KNN method (Ram and Sinha, 2019) to find the nearest trajectory point $t_j^{\prime \prime }$ on T^′′ for each point of the global map $G=\mathit{\{}{g}_{\mathrm{1}},g_{\mathrm{2}},g_{\mathrm{3}},\mathrm{\cdots },g_k,\mathrm{\cdots },g_K\mathit{\}}$. Points of G that share the same nearest point are divided into the same unit $U_j=\mathit{\{}{u}_{j\mathrm{1}},u_{j\mathrm{2}},u_{j\mathrm{3}},\mathrm{\cdots },u_{jm},\mathrm{\cdots },u_{jM}\mathit{\}}$ and marked with the same index j. It is noted that the number M of points of each unit is unequal. Finally, we reorder all points of the global map according to the index j to ensure that the points within the same unit are sequentially adjacent.

The point $t_j^{\prime }$ of the trajectory T^′, which recombines different points of the global map to the same unit U_j, can be used to determine the unit's translation vector A_j. The modified trajectories $T=\mathit{\{}{\stackrel{\mathrm{̃}}{t}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{t}}_{\mathrm{2}},{\stackrel{\mathrm{̃}}{t}}_{\mathrm{3}},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_i,\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_I\mathit{\}}$ (Eq. 6), ${\stackrel{\mathrm{̃}}{T}}^{\prime }=\mathit{\{}{\stackrel{\mathrm{̃}}{t}}_{\mathrm{1}}^{\prime },{\stackrel{\mathrm{̃}}{t}}_{\mathrm{2}}^{\prime },{\stackrel{\mathrm{̃}}{t}}_{\mathrm{3}}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_j^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_J^{\prime }\mathit{\}}$ (Eq. 7), and ${\stackrel{\mathrm{̃}}{T}}^{\prime \prime }=\mathit{\{}{\stackrel{\mathrm{̃}}{t}}_{\mathrm{1}}^{\prime \prime },{\stackrel{\mathrm{̃}}{t}}_{\mathrm{2}}^{\prime \prime },{\stackrel{\mathrm{̃}}{t}}_{\mathrm{3}}^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_j^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_J^{\prime \prime }\mathit{\}}$ (Eq. 8) can be calculated directly by multiplying coordinates of all trajectory points by scaling factors.

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{t}}_i=t_i\cdot \left[\begin{array}{ccc}{f}_{\mathrm{h}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& f_{\mathrm{h}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& f_{\mathrm{e}}\end{array}\right]\text{(7)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{t}}_j^{\prime }=t_j^{\prime }\cdot \left[\begin{array}{ccc}{f}_{\mathrm{h}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& f_{\mathrm{h}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& f_{\mathrm{e}}\end{array}\right]\text{(8)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{t}}_j^{\prime \prime }=t_j^{\prime \prime }\cdot \left[\begin{array}{ccc}{f}_{\mathrm{h}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& f_{\mathrm{h}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& f_{\mathrm{e}}\end{array}\right]\end{array}$$

The translation vector A_j is equal to the coordinate difference of the trajectory point before and after modification (Eq. 9). After moving all the points in the unit with the same translation vector A_j, we attain the corrected unit ${\stackrel{\mathrm{̃}}{U}}_j=\mathit{\{}{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{1}},{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{2}},{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{3}},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jm},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jM}\mathit{\}}$ (Eq. 10) and the adjusted global map $\stackrel{\mathrm{̃}}{G}$, which is the sum of all ${\stackrel{\mathrm{̃}}{U}}_i$ (Eq. 11).

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{A}}_j={\stackrel{\mathrm{̃}}{t}}_j^{\prime }-t_j^{\prime }\text{(10)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{u}}}_{jm}=u_{jm}-{\mathit{A}}_j\text{(11)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{G}=\left\{{\stackrel{\mathrm{̃}}{U}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{U}}_{\mathrm{2}},{\stackrel{\mathrm{̃}}{U}}_{\mathrm{3}},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{U}}_j,\mathrm{\cdots },{\stackrel{\mathrm{̃}}{U}}_J\right\}\end{array}$$

3.2.3 Point cloud smoothing

The modified global point cloud $\stackrel{\mathrm{̃}}{G}=\mathit{\{}{\stackrel{\mathrm{̃}}{g}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{g}}_{\mathrm{2}},u{g}_{\mathrm{3}},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_k,\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_K\mathit{\}}$ is blurry and has a lot of noise because of inevitable oscillations during data collection in the uneven channel. This poses a major obstacle to the observation of internal channel topography and the identification and analysis of channel deposits. Based on the voxel grid (VG) filtering method (Rusu and Cousins, 2011), we propose a new smoothing algorithm, weighted-elastic voxel grid (WEVG), which expands the VG by adding a weight factor and making the voxel size flexible.

The VG filter creates a 3D voxel raster from the input point cloud and then calculates the centroid of each voxel utilizing all points within it. All points within the same voxel are only represented by the centroid. VG filtering can preserve the point cloud's geometric structure (Liu and Zhong, 2014). It is suitable for point cloud maps of the debris flow channel. However, the disadvantages of VG filtering cannot be neglected:

1. The calculation of the centroid is inaccurate without considering point density. In the precise computation of the centroid, density is an independent variable and will affect calculation in an inhomogeneous case.

2. The size of each voxel is the same and fixed so that it is stiff to balance the elimination of outliers and the maintenance of details. The ideal voxel raster should be flexible to partition larger voxels where points are sparse and smaller voxels where points are dense.

The point distribution in the global map is inhomogeneous in this study. The map can be more satisfactory after inserting point density as a weight to the centroid calculation and transforming voxel size adaptive to the distribution of point clouds. We express ρ_k (Eq. 12), the density in point ${\stackrel{\mathrm{̃}}{g}}_k$, as the reciprocal of volume V_k taken up by it, where the volume is defined by applying KNN to search ${\stackrel{\mathrm{̃}}{g}}_k$'s several nearest neighbours and calculating their average distance ${\stackrel{\mathrm{‾}}{r}}_k$.

$$\begin{array}{}\text{(12)}& {\mathit{\rho }}_k=\mathrm{1}/V_k=\mathrm{1}/\left\{(\mathrm{4}/\mathrm{3})\mathit{\pi }{\stackrel{\mathrm{‾}}{r}}_k^{\mathrm{3}}\right\}=\mathrm{3}/\left(\mathrm{4}\mathit{\pi }{\stackrel{\mathrm{‾}}{r}}_k^{\mathrm{3}}\right)\end{array}$$

The voxel is divided based on the principle that each voxel contains the same number 2 N of points, by traversing each point in the global map and searching for 2 N nearest neighbours, recorded as ${\stackrel{\mathrm{̃}}{g}}_{k-N},\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_k,\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_{k+N}$. This voxel division method is extremely fast and can preserve all points with more details while mitigating the influence of the outliers. After calculating the weighted centroid p_k (Eq. 13) of each elastic voxel, we obtain the smoothed global map ${\stackrel{\mathrm{̃}}{G}}^{\prime }=\mathit{\{}{\stackrel{\mathrm{̃}}{g}}_{\mathrm{1}}^{\prime },{\stackrel{\mathrm{̃}}{g}}_{\mathrm{2}}^{\prime },{\stackrel{\mathrm{̃}}{g}}_{\mathrm{3}}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_K^{\prime }\mathit{\}}$ (Eq. 14), which is the aggregation of all centroids, where the WEVG operation is represented by the notation $\widehat{w}$.

$$\begin{array}{}\text{(13)}& {\displaystyle }{p}_k=\left(\sum _{l=k-N}^{k+N}{\mathit{\rho }}_l{\stackrel{\mathrm{̃}}{g}}_l\right)/\left(\sum _{l=k-N}^{k+N}{\mathit{\rho }}_l\right)\text{(14)}& {\displaystyle }\begin{array}{rl}{\stackrel{\mathrm{̃}}{G}}^{\prime }& =\left\{{\stackrel{\mathrm{̃}}{g}}_{\mathrm{1}}^{\prime },{\stackrel{\mathrm{̃}}{g}}_{\mathrm{2}}^{\prime },{\stackrel{\mathrm{̃}}{g}}_{\mathrm{3}}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_k^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{g}}_K^{\prime }\right\}\\ & =\left\{p_{\mathrm{1}},p_{\mathrm{2}},p_{\mathrm{3}},\mathrm{\cdots },p_k,\mathrm{\cdots },p_K\right\}=\widehat{w}\left[\stackrel{\mathrm{̃}}{G}\right]\end{array}\end{array}$$

3.2.4 Cross-section extraction

The cross-sections are vitally important for the development and dynamics of channelized debris flow. The typical cross-sections can also be used to back-calculate the peak flow discharge of debris flow (Whipple, 1997; Berti et al., 1999). Hence, we have developed an algorithm to extract cross-sections of the channel to represent the channel structure more intuitively. Cross-section extraction consists of two parts: cutting and reconstruction.

The channel resembles a meandering curve, so the ideal cross-section cutting plane is the normal plane of each point in the curve. It is feasible to extract the channel cross-sections due to the grouping and sorting of the global map in Sect. 3.2.2. Each single unit ${\stackrel{\mathrm{̃}}{U}}_j$ can be regarded as a cross-section element. Moreover, the smoothed, densified, projected, and corrected trajectory ${\stackrel{\mathrm{̃}}{T}}^{\prime \prime }$ produced in coordinate modification can also reflect the direction of the channel and guide us to determine the position and normal direction of cross-sections. Because the morphology point cloud is sparse, it requires fusing several adjacent elements into a new point cloud ${\stackrel{\mathrm{̃}}{U}}_j^{\prime }=\mathit{\{}{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{1}}^{\prime },{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{2}}^{\prime },{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{3}}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jm}^{\prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jM}^{\prime }\mathit{\}}$. Then ${\stackrel{\mathrm{̃}}{U}}_j^{\prime }$ is reprojected onto one single plane to obtain a 2D cross-section ${\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime }=\mathit{\{}{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{1}}^{\prime \prime },{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{2}}^{\prime \prime },{\stackrel{\mathrm{̃}}{u}}_{j\mathrm{3}}^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jm}^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{u}}_{jM}^{\prime \prime }\mathit{\}}$. When we combine the nearest 2 a units of ${\stackrel{\mathrm{̃}}{U}}_j$ and obtain the new cloud ${\stackrel{\mathrm{̃}}{U}}_j^{\prime }$, we can connect it to the division index trajectory points, ${\stackrel{\mathrm{̃}}{t}}_{j-a}^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_j^{\prime \prime },\mathrm{\cdots },{\stackrel{\mathrm{̃}}{t}}_{j+a}^{\prime \prime }$. The trajectory point ${\stackrel{\mathrm{̃}}{t}}_j^{\prime \prime }$ is located on the projection plane. The normal line of the projection plane is ${\stackrel{\mathrm{̃}}{t}}_{j+a}^{\prime \prime }-{\stackrel{\mathrm{̃}}{t}}_{j-a}^{\prime \prime }$. Then we can determine all points of ${\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime }$ by establishing and solving a parametric equation (Eq. 15), where a new parameter α (Eq. 16) was introduced here.

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\stackrel{\mathrm{̃}}{u}}_{jm}^{\prime \prime }={\stackrel{\mathrm{̃}}{u}}_{jm}^{\prime }-\mathit{\alpha }\left({\stackrel{\mathrm{̃}}{t}}_{j+a}^{\prime \prime }-{\stackrel{\mathrm{̃}}{t}}_{j-a}^{\prime \prime }\right)\text{(16)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }=\left({\stackrel{\mathrm{̃}}{u}}_{jm}^{\prime }-{\stackrel{\mathrm{̃}}{t}}_j^{\prime \prime }\right){\left({\stackrel{\mathrm{̃}}{t}}_{j+a}^{\prime \prime }-{\stackrel{\mathrm{̃}}{t}}_{j-a}^{\prime \prime }\right)}^T/∥{\stackrel{\mathrm{̃}}{t}}_{j+a}^{\prime \prime }-{\stackrel{\mathrm{̃}}{t}}_{j-a}^{\prime \prime }∥\end{array}$$

However, the new cross-section ${\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime }$ is not a curve with a clear boundary, but it is accompanied by a lot of fluctuating noise, which are not qualified for further study. Furthermore, there are many blanks caused by water surface refraction and rock obscuration in the global map, which might make cross-section discontinuous. This is because during data acquisition, the lidar scanning may not cover certain sections of the channel, either because they are unreachable or obscured by the harsh channel environment. These deficiencies impede the analysis of channel structure for subsequent application. Therefore, we smooth cross-section point clouds again with our new smoothing method, WEVG, and then densify it by linear interpolation to gain the final cross-section ${\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime \prime }$ (Eq. 17). Before the linear interpolation, we sort the cross-section in the order of height for wall points and width for ground points.

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime \prime }=\widehat{L}(\mathrm{sort}\left|\widehat{w}\left[{\stackrel{\mathrm{̃}}{U}}_j^{\prime \prime }\right]\right|\end{array}$$

4 Outputs of AscDAMs

The morphological map generated by AscDAMs contains detailed channel morphology data. Section 4.1 explains the selection criteria and the overall characteristics of the optimal maps. Detailed channel morphology including the typical cross-section and basal erosion will be introduced in Sect. 4.2. The delineation of channel deposit distinction, volume estimation, and change monitoring is presented in Sect. 4.3. The map dataset of AscDAMs is on Figshare (Wang and Lu, 2024).

4.1 The global map of the channel

Mapping results calculated with different parameters of downsampling and loop closure will be various. Finding appropriate parameters is a prerequisite for successful and accurate mapping. Among these parameters, the downsampling rate is an uncertainty factor, while the remaining factors either possess predetermined values or exert negligible influence on the mapping accuracy.

After obtaining the trajectory and global map, we projected them onto the DOM of the ZY03 satellite with their elevation data, as shown in Fig. 5, where the relative elevation in the graph is based on the location of the start point of the trajectory (Fig. 5a). The optimal point cloud global map is shown in Fig. 5b. Detailed channel morphology represented by different typical cross-sections (TCSs) (Fig. 5b) is selected for further demonstration.

Figure 5 Overview of the AscDAMs map, projected on the DOM of ZY03 satellite: (a) the trajectory and (b) the global map. TCS A–D are typical cross-sections of the channel.

The uneven trajectory (Fig. 5a) shows the ruggedness of channel interior and the difficulty of employing original LIOSAM to such a harsh environment. It is clear to see that the trajectory we calculated and channel area on the DOM are highly matching. The global map (Fig. 5b) presents the same good fitting. The structure, outline, and confluences of the channel point cloud obviously match well with the satellite image, although there is a slight offset at the upstream end. We can also see that the height of the cliff captured in the map is tens of metres (Fig. 5b). The maximum elevation difference of the channel bed within the trajectory range is about 160 m. Two drifts on the trajectory seem incongruent with the DOM and will be discussed in Sect. 5.1.3.

The DEM with a horizontal resolution and elevation accuracy of 0.1 m was derived from the global map and is shown in Fig. 6a. It should be noted that the DEM generated from AscDAMs was based on the filtered point cloud global map. Because of the existence of basal erosion and an overhanging cliff, there could be several elevation values at the same raster. Hence, the redundant elevation values should be removed to avoid unnecessary errors, although this filtering process will lead to a decrease in channel morphology information. The DEM obtained by AscDAMs can be integrated with a DEM generated by satellites and/or UAVs, which can be used in numerical simulations of debris flows.

Figure 6 The channel DEM and coloured map: (a) the AscDAMs DEM of the channel, (b) withered grass, (c) deposit and fresh plant, and (d) bedrock in shadow.

The coloured information is also integrated with the global map by running RGB–SLAM, enhancing the comprehension of the channel interior information. The channel map is predominantly grey-brown after rendering because the data were collected in February when most of the vegetation did not sprout. In the map, unconsolidated sediment is grey, while withered grass is brownish (Fig. 6b). The long-standing stable deposits are brown (Fig. 6c). There are many places in the channel where the colour is dark green. These sections are usually full of fresh plants (Fig. 6c), while bedrock in the shadow might exhibit the same colour (Fig. 6d). The coloured topographic map can be utilized for determining the deposit stability and vegetation restoration.

4.2 Detailed channel morphology

4.2.1 Typical cross-section

Although the channel morphology varies throughout its entirety, the geomorphic features within the channel can be reflected by different typical cross-sections as shown in Fig. 7, such as artificial facilities (TCS A1), rocks (TCS A2), and pools (TCS A3). Human-made facilities like a small dam are the easiest to identify due to their regular shape, while a natural landscape needs more attentiveness. By cutting out and densifying cross-sections, the various shapes of the channel are more intuitive. Typically, the cross-sections of the channel are V-shaped or U-shaped, featuring a narrow bottom surface, wide top, and steep side walls.

Figure 7 Inside view of TCS A1, A2, and A3 in the smoothed map (a, d, g), camera images (b, e, h), and cross-sections before and after reconstruction (c, f, i).

At the channel bed, the ground is uneven, with deposit materials having a wide grain size distribution from silt or clay to cobbles and huge rocks. Since the water bodies will reflect and refract the lidar signal, there are some missing parts in the point cloud. For distinguishing blanks, narrow, continuous, and low-lying areas on the ground are mostly streams, while large horizontal gaps are mostly pools. The other types of blanks are mostly due to missing scans. On both sides of the channel, the sleek channel boundaries correspond to clear and smooth point clusters, while vegetation and loose material correspond to fuzzy and rough clouds. The distinct and detailed morphology at the bottom of the channel, including the channel bed and side walls, which is significant to the development of debris flows, cannot be obtained from existing satellite images.

Figure 8 The erosion pattern in TCS B1 (a–c) and B2 (d–f) with a large trajectory drift, as shown in Fig. 5a.

Figure 9 The deposit in TCS C and volume estimation.

4.2.2 Erosion pattern

Debris flow erosion consists of two parts: bed entrainment and bank erosion (Stock and Dietrich, 2006; Berger et al., 2011a). Bank erosion is important for channelized debris flow, especially when debris flow changes its direction. Satellite images could not provide sufficient channel interior information for the study of erosion. However, the present study successfully addressed this bottleneck technique question using the global map and extracted cross-sections. As an illustration, AscDAMs results distinctly reveal two areas of bank erosion, as shown in Fig. 8. The cliffs lean towards the channel interior because their middle and lower parts have been eroded heavily. These erosions can be caused by the impingement of debris flows and/or other friction. From the perspective of the DOM, the shadow on these two segments is heavy. It indicates that the cliffs are towering and steep so that the environment inside the channel is difficult to detect. Thus, the AscDAMs holds a unique advantage in detecting lateral bank erosion and providing supplementary information that cannot be observed from satellite images and UAVs.

4.3 Channel deposit distinction and monitoring

4.3.1 Channel deposit distinction and volume estimation

Deposit in the channel can be identified directly by the global map and cross-sections. At the channel bed, there are scattered stones in most places, which are either uneven on the surface, such as TCS A2 (Fig. 7), or inclined to accumulate at the foot of the sidewall, such as TCS B1 (Fig. 8). Besides, individual large deposits can also be clearly recorded by the global map. For instance, at TCS C (Fig. 9), a body of deposit can be identified of about 7 m high, occupying about half of the channel bed width, next to the north channel wall. The erosion pattern on this deposit clearly identified that the part outside of the bend was severely eroded. Furthermore, the volume of deposits can be estimated with the coordinate data of the global map and cross-sections. The deposit in TCS C, for example, can be approximated as a combination of regular polyhedrons, as shown in Fig. 10a. Assuming the bedrock surface behind the deposit is flat and regular, we can divide the deposit body of TCS C into a quadrilateral prism with a trapezoidal base and a triangular prism with a trapezoidal side. The volumes of the two prisms are 661 and 504 m³ respectively. Hence, the total volume is estimated to be 1165 m³. Using the same method, the deposit volume in TCS A2 (Fig. 10b) and TCS B1 (Fig. 10c) is evaluated to be 540 and 595 m³ respectively. Note that this is a rough estimation as the deposit boundary inside is not known. If the channel morphology is detected by AscDAMs before and after the debris flow event, the accuracy of estimating volume change could be significantly improved.

Figure 10 Deposit volume estimation in TCS C, TCS A2, and TCS B1, approximated as regular polyhedrons.

Figure 11 The deposit changes after collapsing in TCS D1.

A proper estimation of possible source material of debris flow is vital to estimate the volume and destructive power of debris flow, which are significant to hazard mitigation. However, the deposits are difficult to identify and evaluate from the DOM or the DEM of satellites, especially in narrow alpine channels. Satisfactorily, it can be identified intuitively and analysed quantitatively on the global map now, which demonstrates the technological advantage of AscDAMs.

4.3.2 Deposit change monitoring

Regular detection of channel morphology by AscDAMs makes it possible to study the long-term spatial and temporal evolution of channel deposits. In the current study, two field investigations are carried out in February and November 2023 which span an entire rainy season. A new dataset was obtained by AscDAMs after a whole rainy season. By comparing the two different data, it is possible to monitor the migration of loose channel deposits during the rainy season. The results of AscDAMs show that the deposits inside the channel almost have no significant change during this period, except for a collapse of the deposit at TCS D1, as shown in Fig. 11. The deposit was like a triangular prism tightly against the channel wall last time. After collapsing, it resembles a lying pyramid, with negligible change in volume. The quantitative analysis of changes in channel deposits provides valuable insights for studying debris flow risk in the study area. This aspect holds significant importance for hazard mitigation efforts.

Figure 12 Point clouds of TCS D1 (a, b), D2 (c, d), and D3 (e, f) before and after smoothing.

5 Discussion

5.1 The quality of smoothing and extraction

The quality of the point cloud smoothing can be exemplified by TCS D1, D2, and D3 (Fig. 12). The position of TCS D1 is a confluence point of the tributary to the main channel of the Chutou gully (Fig. 5b). A small alluvial fan was formed at the outlet of the tributary. The surface of the alluvial fan and rocks in the smoothed images is clearer than that in the initial map. The smoothed point cloud of an artificial wire is clearer and more distinguishable. TCS D2 is another intersection. A 300 m³ deposit's outline in front of a tributary outlet is more distinguishable after smoothing. There is a pool at TCS D3. The boundary of the pool is more distinct after smoothing, and there is significantly less noise on the adjacent rock wall. Compared with the point cloud images before smoothing, the visual effect after smoothing is significantly improved with clearer material boundaries and a lighter blurriness. It makes observing the landforms inside the channel much more convenient and accurate.

The quality of cross-section extraction is also satisfactory. In cross-section extraction, both images before and after reconstruction are preserved (Figs. 7–9). Point clouds of cross-sections before reconstruction are sparse and cluttered, while the reconstructed ones are dense and ordered. Even if a reconstructed cross-section has some offset on occasion (Fig. 9c), it does not affect its description of the channel profile. For comparison, we created cross-sections of TCS A1 and TCS B1 with the data obtained by a laser range finder (LRF) from field investigation and the DEM from ZY03 stereo images, as shown in Fig. 13. Compared with the LRF and ZY03's DEM data, the reconstructed cross-sections and the AscDAMs DEM are more accurate, with more channel structure information. Both ZY03's DEM and the LRF data lose channel details. Moreover, the errors of ZY03's DEM are apparently large and cannot reflect lateral erosion of the channel inner wall and erosion of deposits. In general, the AscDAMs results are superior to the commonly used data from existing methods. Their accuracy is much higher than that of the prevailing satellite-derived DEM. Such high-quality cross-sections and DEMs can be further used for simulating and even predicting debris flow. This precision is crucial for a comprehensive understanding of debris flow mechanisms.

Figure 13 Cross-sections from ZY03's DEM, the LRF, the AscDAMs DEM, and global map before and after reconstruction: (a) TCS A1 and (b) TCS B1.

5.2 The accuracy of AscDAMs

5.2.1 The elevation error of AscDAMs

We analyse the elevation error by comparing the height data of the initial trajectory obtained by LIOSAM, the modified trajectory optimized by AscDAMs, barometer data, and the DEM of ZY03 (Fig. 14). Taking the barometer result as a benchmark, the average bias and root mean square errors (RMSEs) of the LIOSAM, AscDAMs, and DEM trajectories are estimated through polynomial fitting of the height data (Table 3). Evidently, the DEM of ZY03 exhibits significant errors within the research area, featuring abnormal sudden rises and drops of elevation that deviate from the field investigations. The fluctuation of elevations of ZY03 can reach 100 m (Fig. 14). It is found that the AscDAMs trajectory has the most precise elevation data with an average bias of −0.48 m and a root mean square error of 8.65 m. This affirms the effectiveness of the proposed deviation correction algorithm. The satellite-derived elevation, namely the DEM of ZY03, is incorrect, with many unreal fluctuations along the channel (Fig. 14), which is due to the inability of the satellite to detect the channel bed where it is narrow and in shadow (Cao et al., 2021).

Figure 14 The height data of the initial trajectory obtained by LIOSAM, the barometer, the DEM of ZY03 on 14 January 2021, and the modified trajectory optimized by AscDAMs.

Table 3 The average elevation bias and RMSE of the LIOSAM, AscDAMs, and DEM trajectories.

5.2.2 The horizontal error of AscDAMs

The bias, displacement error, and distance error between the end point of the trajectories and the referred end point are calculated respectively based on the DOM (Fig. 15), as shown in Table 4. The reference displacement, which is the straight distance between the referred start and end points in the DOM, is 1090 m. The reference distance, which is the estimated length of the channel from the referred start point and referred end point, is 1190 m. The horizontal bias of AscDAMs is 28.66 m, resulting in a displacement error and a distance error of only 2.63 % and 2.43 % respectively. As a comparison, the initial maps of LIOSAM have a large offset of 182.68 m. This shows that the proposed deviation correction algorithm in AscDAMs can significantly improve the horizontal accuracy from the original LIOSAM result.

Figure 15 The trajectories of AscDAMs after deviation correction and LIOSAM (i.e. the initial trajectory).

Table 4 Horizontal bias of AscDAMs and LIOSAM.

5.2.3 Drifts and bank erosion

As mentioned above, there are two irregular deviations of trajectory: Drift 1 and Drift 2 (Fig. 5a). However, these two drifts are not introduced by the accumulation of systematic error. By comparing the locations of drifts and selected bank erosion areas, it is easy to find that the drifts correspond to the bank erosion areas. The obscuration by overhanging cliffs prevents the DOM from capturing accurate channel morphology. Consequently, the real channel morphology could not accurately be depicted on the DOM of ZY03, contributing to the confusion of “drifts”. This shows the superiority of AscDAMs to capture details of channel interior structure over satellite images.

5.3 The universality and usability of AscDAMs

AscDAMs succeeds in promoting SLAM to complex gully environments to obtain high-resolution topographic maps with full characteristics of debris flow channels. In the previous study, SLAM has been tested in a small-scale hillslope gully with gentle terrain (Kinsey-Henderson et al., 2021) and broad hillsides with stable operating conditions for SLAM (Marotta et al., 2021). In the current study, the research area was selected as the most challenging environment for AscDAMs. The DOM-aided deviation correction algorithm effectively minimized the error and greatly enhanced the accuracy of the SLAM results. The point cloud smoothing algorithm mitigated the effects of sensors swaying and rocking, leading to a substantial improvement in the quality of the final point cloud map. By utilizing the cross-section extraction algorithm, the channel morphology could be quantitatively assessed. The location and volume of channel deposits can be precisely quantified. The change of each deposit could also be accurately detected by comparing two different AscDAMs maps. The successful implementation of AscDAMs in such a challenging environment implies that, for other channels in the Wenchuan earthquake region or less demanding scenarios with more typical features for the SLAM algorithms to compute, favourable results can also be achieved.

AscDAMs is a multi-sensor fusion system. Compared to existing channel detection technologies, AscDAMs is easy, economical, and effective. It can be assembled with only three core components: lidar, the IMU, and a minicomputer. Complementary to the current field investigation methods, our system can record the channel morphology quantitatively and automatically while travelling along the channel without manual supervision or intervention. AscDAMs does not need complex or elaborate route planning in advance. Furthermore, without any special requirement of hardware, an ordinary computer is sufficient for this calculation, and the computing time is equivalent to the data collection duration. If a higher accuracy is pursued, the global map can be further optimized with input from magnetic sensors, altimeters, and other sensors. AscDAMs can utilize different sensors according to specific environments, while the underlying algorithm's logic and processing flow remain unchanged.

5.4 The limitations of AscDAMs

Although AscDAMs offers advantages in terms of accuracy, novelty, and efficiency, it does exhibit certain limitations to some extent.

1. The final global map will become sparse visually after smoothing for point cloud resampling (Liu and Zhong, 2014). This is the common defect of point cloud noise filtering and smoothing algorithms.

2. The colour contrast of the coloured map is insufficient. Vegetation changes with season. In this study, we only tested AscDAMs before spring when plants had not yet recovered. More tests can be implemented in different seasons. In addition, the complexity of the light affects the colour recognition of the camera.

3. The equipment is applicable for areas within the reach of human resources. In some cases, channels might be blocked for walking. AscDAMs can be combined with more topography-adaptive carriers in the future.

4. The systematic error inherent in the SLAM algorithm has not been entirely overcome. Although the proposed deviation correction algorithm significantly mitigates the systematic error of the SLAM algorithm, it remains impossible to completely eliminate the influence of this error.

6 Conclusion and perspective

Obtaining the high-resolution, accurate topographic channel map is the common key challenge for channelized debris flow research. At present, widely used satellite images, UAV-based mapping, and other existing technologies cannot satisfy the requirements of accuracy and efficiency in observing channel interior conditions in long and deep mountainous gullies. SLAM is an emerging 3D mapping technology and has been applied across different platforms for numerous scenarios of topographic mapping. However, state-of-the-art SLAM mapping results contain large drift and abundant noise induced by the extremely rugged, long and deep channel environment. Aiming to solve these problems, we proposed AscDAMs with a set of new algorithms, including deviation correction, point cloud smoothing, and cross-section reconstruction to process the original SLAM results. In addition, a map colouring algorithm is developed to supplement more information and the map. A frequent debris flow gully named Chutou gully in the Wenchuan earthquake region was selected as the research area. AscDAMs was successfully implemented in extremely harsh environments, resulting in the high-resolution full-character morphological mapping of debris flow gullies. Compared to existing channel detection technologies, AscDAMs offers the following benefits:

1. Improved accuracy. The proposed deviation correction and point cloud smoothing algorithms significantly enhance the accuracy of mapping results.

2. Cross-section extraction. The cross-section extraction algorithm enables the full characterization of debris flow channel cross-sections, facilitating the study of critical channel cross-sections in terms of debris flow development, dynamics, and erosion.

3. Comprehensive 3D mapping. The 3D mapping with adequately detailed information is sufficient to quantitatively assess the position and spatial distribution of channel deposits. With reasonable simplification, it is also possible to estimate the deposition volume, which is vital for risk assessment and management.

4. Morphological monitoring. Periodic resurveys of the channel with AscDAMs enable the monitoring of gully morphology changes, such as the downward movement of loose slope material and sediment transport within the channel, from a quantitative standpoint.

5. Vegetation recovery analysis. The additional colour information captured can be utilized to study vegetation recovery inside the channel.

As a crucial supplement to existing channel morphology detection methods, AscDAMs works well in the complex channel environments. It provides the important but currently absent channel interior details, which is promising to promote deep understanding of debris flow mechanisms and post-seismic long-term evolution. It can also support precise hazard–risk assessment and mitigation, although it can be further improved in systemic error correction.