A new explainable DTM generation algorithm with airborne LIDAR data: grounds are smoothly connected eventually

3-dimensional coordinates of the earth’s surface are usually collected as point cloud data, and the point cloud is assumed to be proper enough to model DTM. However, the point cloud is still a set of observations of real-world entities, so it necessarily has a limitation in perfectly representing the world. In particular, sensors most widely used for DTM generation, including LiDAR, have limitations in that their outputs are irregularly spaced 3-dimensional coordinates and their local point densities are inevitably varying. Even if the same laser pulse rate has been used during the flight mission, the point spacing will inevitably be different as it is affected by many factors such as flight configuration, flight condition, and objects on the ground (Habib et al., 2011; Morsdorf et al., 2008; Yu et al., 2004). This is the main reason why most algorithms based on point clouds or TIN representations require lots of parameters to tune and are hard to generalize in universal scenarios.

In contrast, the image grid representation provides regularly spaced data. Specifically, when point clouds are the observation of the earth’s surface from airborne laser scanning, the image grid representation of point clouds is DSM. To generate DSM, users need to determine the grid size for rasterization, and the grid size often becomes the resolution of DTM. Depending on the grid size, multiple point clouds can share the same grid and some grids might not have any points. The down-sampling, which generates a large coarse grid DSM, often alongs with the rasterization to prevent void grids (Hyyppa et al., 2001; Maltezos et al., 2018; Oh et al., 2022). Here, this type of rasterization with down-sampling is referred to as a “coarse rasterization”. The coarse rasterization can prevent void grids but results in data loss. To alleviate the data loss problem, this study adopts a “fine rasterization” that projects point clouds into a finely and regularly spaced image grid.

Figure 1 provides a graphical illustration comparing coarse rasterization and fine rasterization. To project 7 observation points into an image grid, the coarse rasterization uses an image of 2 by 2 grids while the fine rasterization uses an image of 3 by 3 grids. The coarse rasterization does not have void grids while the fine rasterization has void grids initially. Void grids of the up-sampled image grid in the fine rasterization are filled up with the nearest neighbor interpolation. Fine rasterization can preserve more observation points with smaller displacement (registration) error than coarse rasterization, and thus it can help to generate precise, high-resolution DTM. However, the image grid representation either with coarse or fine rasterization may still result in data loss if some points coexist on the same grid. When multiple points occupy the same pixel, our DTM generation algorithm uses the lowest elevation point (the last return of LiDAR points) for the DSM value as the ground typically be the lowest elevation among neighboring points.

A Sobel operator is a widely used kernel, particularly for edge detection in image processing applications as it essentially computes the gradient of the image intensity (Abdou and Pratt, 1979). Typically, a Sobel operator convolves two 3 by 3 kernels with the original image where kernels calculate derivatives of horizontal and vertical directions, respectively. Then, the Euclidean norm of two derivatives can describe the gradient of each pixel of the image. When the image is DSM, the gradient can be used to approximate the slope of the surface (Gelbman and Papo, 1984). Therefore, DSM can be transformed into a slope map that describes the slope of the topography. Based on the slope map, we delineate a break-line map that shows the line where the topography is steeper than a certain level of slope. Figure 2 (a-c) illustrates an exemplary procedure of break-line mapping with the Sobel operator.

A common assumption for DTM generation is that ground is generally smooth while non-ground objects have protruding shapes. In addition to this conventional assumption, we add the assumption that any non-ground object is surrounded by steep slopes while grounds are smoothly connected eventually. Thus, the proposed algorithm filters out any area surrounded by more than a certain degree of slope (i.e., the slope threshold). This assumption is reasonable and robust as hills, high-relief terrains, cliffs, mountain ranges, valleys, and overpasses are eventually connected to smooth surfaces in most cases; while non-ground objects like buildings and trees are enclosed by steep slopes or break-lines. In addition, since the break-line can span infinitely, the algorithm is not confined to the local neighborhood but can consider the global neighborhood and can classify non-ground objects regardless of their sizes and shapes. As a result, a non-ground object in our algorithm is clearly defined as a set of break-lines and an area surrounded by break-lines. Lastly, pixels classified as non-ground objects are masked and linearly interpolated with neighboring ground elevations to produce a seamless rasterized DTM.

The essential parameter of the proposed algorithm is the slope threshold parameter which is in charge of delineating break-lines. We claim that the parameter clearly possesses a physically meaningful value, and it can be easily tuned based on the user’s objective and topographical characteristics. We set the slope threshold of 45 degrees as default because it is robust enough to produce reliable DTM in most topographies. The impact and suggestions for parameter selection are discussed in Section 3.2.1.

A few additional considerations were put into the algorithm to increase its scalability and facilitate its practical usage. First, some non-ground objects can lie on the edge of any given DSM layer. In this case, the non-ground object cannot be classified as a non-ground object because it is not fully surrounded by the break-line but is partially opened due to the limited data extent (the extent of given DSM layer). Therefore, the algorithm set the edge of the DSM as a break-line initially. Note that this action will enclose the ground that was originally not enclosed with break-lines. To prevent this issue, a decision based on the area was made to determine whether an enclosed area is ground or non-ground: First, if the enclosed area is smaller than a low limit (A1), the enclosed area is determined as a non-ground. Second, if the enclosed area is larger than a high limit (A2), it is determined as ground. Third, for the enclosed area between a low limit (A1) and a high limit (A2), a metric called the rectangularity that describes the ratio of the enclosed area to its minimum bounding rectangle area was considered. If the rectangularity is larger than a certain value (R), the enclosed area is classified as a non-ground, otherwise, it is classified as ground. This decision is based on the assumption that relatively large non-ground objects are mostly large buildings with rectangular shapes. Also, as all areas are eventually bounded by the size of the data extent, A2 is necessary. We set A1, A2, and R as 40,000$m^{2}$, 100,000$m^{2}$, and 50% as a default. This is because objects larger than 40,000$m^{2}$ are rectangular-shaped buildings in most cases. We found these default values hardly produce artifacts. Figure 2 (d-e) displays an exemplary procedure of non-ground filtering, and Figure 2 (f) shows the final DTM. Figure 3 illustrates a ground and non-ground classification rule of the proposed DTM generation method.

As the main task of DTM generation is typically considered as a classification of ground and non-ground, a group of algorithms called the “ground filtering algorithm” has received lots of attention instead of DTM generation (Meng et al., 2010). Also, studies often evaluate the performance of DTM generation based on several binary classification metrics (Sithole and Vosselman, 2004; Mongus and Žalik, 2013; Hu et al., 2015). However, the ground filtering algorithm is usually limited in its purpose for classifying ground and non-ground by its definition and is not for mapping a full extent of a digital map that include both grounds and water bodies. Also, in general, external sources for water mapping are readily available due to lots of accumulated remotely sensed imagery (Huang et al., 2018). Perhaps these are the reasons why most DTM generation algorithms have ignored water body mapping (Hu and Yuan, 2016; Gevaert et al., 2018; Amirkolaee et al., 2022) even if subsequent analyses based on DTM often require a water map. However, the use of external data sources can lead to errors due to registration issues or mismatches in spatial and temporal resolution with the LiDAR data. Also, obscuration from clouds can be a problem when timely mapping is needed. In addition, a water map itself can be helpful to prevent errors in DTM mapping (Susaki, 2012).

Therefore, we include a function that can extract water bodies and their elevations as well into our open-source workflow of DTM generation. In the proposed algorithm, water pixels are identified based on the assumption that the point density over water bodies is much lower than in non-water areas as water bodies hardly reflect laser points. With the finely rasterized DSM before the interpolation, the average point density (P) of a given scan is calculated by dividing the number of non-void grids by the number of total grids. Then, the number of non-void grids in a certain size of a sliding window (N pixels) will have a binomial distribution B(N, P). Specifically, B(N, P/2) was used to compensate for the imbalance of point density due to scanning overlap and to avoid overly detecting water. Based on the binomial distribution, lower confidence bound was used for the decision boundary of water classification. We used a window of 9 by 9 and a confidence level of 4 as a default. Lastly, the elevation of water bodies was selected by taking the 10th percentile of elevations among each water segment to prevent outliers. These parameters were set empirically and found to be robust in diverse topographic airborne laser scanning, but it is worth noting that elevation values cannot guarantee the true elevation as observations of water bodies contain lots of noise. This is because a LiDAR for topographic mapping commonly uses a near-infrared laser, which is absorbed by water and cannot reflect the laser point. Moreover, the elevation of the water bodies is dynamic in nature due to the water cycle. A more detailed description and impacts of water-related parameters are provided in Section 3.2.3.

The proposed DTM generation algorithm can convert a points cloud of airborne laser scanning to a rasterized DTM. The algorithm starts by generating the finely rasterized DSM to keep original points and transform the data regularly gridded. Based on the assumption that all non-ground objects are enclosed by a certain level of a steep slope, the algorithm delineates a break-line map with a Sobel operator and classifies non-ground objects based on the rectangularity and the size of the enclosed area. Finally, water mapping is performed considering the point density. Our method performs in an end-to-end manner and is easy to use as the meanings of parameters are very straightforward. Also, its results and errors are explainable and predictable in general, which can greatly reduce uncertainties in the resulting DTM.

This subsection provides the comparison results among OUR, CSF, and LAS. We excerpted four tiles that show distinctive differences among different methods and that can provide helpful information to potential users.

Figure 5 shows aerial RGB images and three different DTMs from different methods. The RGB images are from the U.S. Department of Agriculture’s (USDA) National Agriculture Imagery Program (NAIP)’s orthoimagery. The rank denoted with RGB images indicates the order of the highest MAE values out of 81 tiles. The elevation ranges of OUR’s DTM were provided for reference. The RMSE and MAE values of other DTMs calculated by comparing to OUR’s DTM were also provided.

Figure 5 (a) displays an urban area along the river. CSF and LAS were not able to generate proper DTM for the large building. CSF regarded the large building as the ground while LAS produced a hole (0 value) as same as the river. On the other hand, OUR filtered out the building as a non-ground and interpolated it with nearby ground elevations. Another unique difference can be found in the bridge. As the bridge is connected to the ground without a discrete elevation change, OUR classified the bridge into a ground category. However, neither CSF nor LAS considered the bridge as a ground object.

Figure 5 (b) shows the overpass structures. As the overpass is connected to the ground smoothly, OUR classified the overpass as ground. However, both CSF and LAS considered the overpass as non-ground. This is because CSF simulates a cloth covering on the upside-down DSM. The TIN-based method, LAS, also produced a similar result. Another difference is that OUR regarded the pile of soil as the ground while CSF and LAS removed it from the ground category.

Figure 5 (c) illustrates the disadvantage of OUR as it overly smoothed deep valleys (i.e., West Lafayette Parks Maintenance). The reason why the small, narrow valleys were smoothed is that some narrow valleys were enclosed by steep terrains.

Figure 5 (d) shows the advantage of OUR in that our method was able to filter out a large building as a non-ground and produced a reliable DTM while a significant portion of the large building remains in both DTMs of CSF and LAS.

In summary, OUR produced more reliable DTMs compared to CSF and LAS, particularly for urban areas with large buildings. Also, OUR has a unique characteristic that can map bridges and overpasses as a category of the ground. It does not mean OUR falsely classified bridges and overpasses just because bridges and overpasses are human-made structures. Although it is not a natural terrain, it is an artificial terrain like roads and plays a role more like a ground. Moreover, the terrain under the bridge and overpasses are actually unknown. In addition, note that both CSF and LAS partially mapped overpasses as ground as shown in Figure 5 (b). CSF and LAS simply interpolated unmeasured grounds after removing the layer on top. Also, considering CSF and LAS classified some parts of the overpass as ground classes, we can argue that CSF and LAS were not consistent in their decision rules for overpasses, and their results are difficult to predict and explain. Rather than debating whether bridges and overpasses are ground or not, it is worth noting that each DTM definition has its own merits. DTM that regarding bridges and overpasses as the terrain can be beneficial in some applications where overpasses and bridges should be considered as ground such as building extraction (Song and Jung, 2022). Lastly, OUR method with the default parameter overly smoothed some steep, narrow areas as shown in Figure 5 (c). However, it can be prevented by tuning the slope threshold. The impact of the slope threshold is detailed in Section 3.2.

Since DTM generation involves the binary classification between non-ground and ground, our method shares the common limitation of binary classification that has a trade-off between omission and commission errors. To be specific, as the slope threshold increases, the DTM of the steep area retains sharper terrain reliefs, but some non-ground objects can remain as ground. Conversely, if the slope threshold decreases, the DTM is getting more smoothed while it could prevent the presence of non-ground objects in the resulting DTM. This limitation commonly exists in DTM generation algorithms and the trade-off can be controlled by several parameter tunings in some algorithms (Chen et al., 2017; Liu, 2008; Meng et al., 2010). For example, CSF can tune the rigidness of cloth (Zhang et al., 2016). LAS can tune the parameter for the maximal standard deviation for planar patches of TIN (Axelsson, 2000). It is worth mentioning that the parameter setting of our method is very straightforward and intuitive, and thus, the outcome according to the parameter is easily predictable, compared to other methods. This advantage enables users to find the proper parameters for their objectives and study areas and can significantly reduce uncertainties in the resulting DTM.

Although the proposed algorithm was confirmed to be robust to generate DTM of diverse topography, our algorithm can suffer when only the physical shape of the target is not enough to identify whether it is a ground or non-ground. This limitation is common in most DTM generation algorithms. For example, there is no way to distinguish between dome-shaped buildings and the same shape of rocky terrain unless a sophisticated semantic understanding is possible. Adding more delicate decision rules might enable the algorithm to distinguish some rare but difficult cases, but it will compromise the algorithm’s generalization capability and computational efficiency. Deep learning-based algorithms that can make a decision based on the learned distribution may work better, particularly when sophisticated semantic understanding is needed. However, it will necessarily entail errors from the distribution shift where the distribution of the target area is different from that of the trained area (Moreno-Torres et al., 2012; Tuia et al., 2016). Also, another limitation of deep learning-based algorithms is their decisions are necessarily limited by their input size (Amirkolaee et al., 2022; Gevaert et al., 2018; Hu and Yuan, 2016). The typical input size of deep learning-based models for semantic segmentation is 256 by 256. Therefore, the model should decide whether the target is ground or non-ground based on the limited area of 256-meter by 256-meter if the resolution of DSM is 1-meter. However, if the input is composed of only the center of a large flat building, the deep model will be likely to fail in its mapping. Errors in large object detection frequently occur in deep learning-based methods (Song and Jung, 2022; Ji et al., 2018). On the other hand, our proposed algorithm can consider the entire data extent for the decision. In other words, inputs do not need to be tiled like in deep learning. More importantly, due to uncertainties of the decision rule of deep learning methods, the so-called “black box”, predicting and explaining the results would be difficult, which could jeopardize the credibility of the subsequent analyses. Of course, our algorithm is not without errors. However, the magnitude and influence of the error can be better estimated than other algorithms as parameter tunings are very intuitive and the result of the algorithm can be easily explained.

Another limitation to be noted is the uncertainties in water elevation and water mapping. Our DTM mapping workflow includes a feature for water body mapping to alleviate elevation errors near water bodies and to assist subsequent studies. Having the water mapping in our workflow is advantageous as users can expect a DTM having all elevations in both terrain and water, and it can replace the post-processing of water mapping that requires another external source of data. However, due to the low reflectance of a near-infrared laser to water, observations include lots of noise, resulting in uncertainties of water elevation. In fact, even if the MAE was calculated after masking the water area when performing the tile comparison in Figure 5, tiles containing large water bodies recorded the highest MAE (18 out of the top 20 contain water bodies either river or lake). This suggests that the water mask was not perfect and that the CSF and LAS also had errors near water bodies. Particularly, we witnessed water with fast flowing streams and lots of floating often has high point density, resulting in being omitted from the water mask. There have been previous studies that map water with airborne LiDAR data by using a supervised classification (Brzank et al., 2008; Smeeckaert et al., 2013) or LiDAR signal intensity (Höfle et al., 2009). Although they require either training procedures (Brzank et al., 2008; Smeeckaert et al., 2013) or a radiometric correction of intensity (Höfle et al., 2009), they could be an alternative way of mapping water bodies. Future studies for more accurate and scalable water mapping with airborne LiDAR data are needed.