Algorithm of Multiple Filter to Extract DSM from LiDAR Data

Masaomi Okagawa

This paper describes the processing method of the LiDAR data by utilizing the feature characteristic derived from height information. The LiDAR data is classified by the neighborhood height effect using clustering process of statistics. It is recognized using the control point made by each cluster whether the cluster is surface objects or non-surface objects (e.g., trees, buildings, cars etc.). On the other hand, the filtering process that added pulse return information is carried out in forest areas and then LiDAR data is classified by the areal characteristics using the functions of ArcView GIS 3.2 and its extensions.


1. Introduction

The use of Light Detection and Ranging (LiDAR) data for terrain models and topographic mapping has gained wide attention in the recent years as it contains the height data on the earth surface. Global Positioning System (GPS) and Inertial Measurement Unit (IMU) acquire the position and the attitude of the airplane, and the laser measurement unit acquires the distance to the earth surface. After that the point coordinate X, Y, Z is obtained by analyzing the position, attitude, distance and projection transformation. Since this data includes many Non Surface Objects such as buildings, trees and cars, the classification is needed in various application fields. This classification is called the filtering process of LiDAR data.

The filtering process, based on the edge extraction method, has been developed for the newly applied technologies 1). However, the area should be divided in detail to set up the proper parameters of the edge extraction especially in Japan where the urban areas have very dense buildings and steep undulating topography exists in forest areas. This process is mainly done manually and it is very time and cost consuming process. Therefore the introduction of a new method to speed up and to automate this process is necessary. The problem of traditional technology is also researched. The interpolation is needed to apply image processing because LiDAR data is a digital point data. This error (Interpolation Error) has an effect on the result. Therefore, the development of the technology that processes raw LiDAR data itself becomes important.

The purpose of this research is to extract DSM (Digital Surface Model) from LiDAR data. Sometimes the raw data itself is called DSM and sometimes not. In this paper, DSM is defined as following to eliminate complications.

In other words, DSM is the model that includes existing terrain and non-surface objects (e.g. trees, buildings, cars etc.).

This research proposes a method that processes the raw data itself by using the concept of the clustering of the statistics and multivariate analysis 2). Clustering is the concept that is applied as a technique of pattern recognition in the field of image processing 3). The goal of this development is to computerize to identify the surface of the earth and/or non-surface objects on the 3D profile of the raw LiDAR data in the same way as a human user does it by visual interpretation. However, the analysis or the recognition that is judged by the human interpretation is faced with the following problems 4).

While solving such problems, the algorithm of DSM generation is discussed. Various software systems were used including ArcView GIS 3.2 and its extensions for the data processing.

2. Algorithm

The filtering algorithm for the DSM generation as defined in Section 1 is discussed here. The property of data is explained before the main subject. The input LiDAR raw data has rectangular coordinates X, Y, Z and return information as the attributes. The pulse that has returned first is called the first pulse, the second pulse, and the third pulse in order, and the pulse that has returned last is called the last pulse. The output is the DSM as defined in the introduction part. The return information becomes more effective in forest areas. Estimation of tree and surface height becomes possible by return information. The laser pulse does not penetrate the branch and leaf but reaches to the ground as the sunbeams streaming through the leaves, because the diameter of the pulse is 20`30 p on the ground. However, the first pulse is not necessarily a reflection on the top of the tree and the last pulse is not a reflection from the ground. Therefore the proper analysis should be processed to utilize return information. This issue will be mentioned in Section 2.4.

2.1 Cluster Analysis

The classification of raw LiDAR data applying the concept of the Cluster Analysis that is used with Statistics and Multivariate Analysis is discussed here. The Cluster Analysis is a method that classifies into some groups that have the resembling characteristics. First of all, raw LiDAR data needs to be transformed into the point in the Characteristic Space. The Characteristic Space is the lower dimensional space than the original pattern that expresses its property accurately. The sufficient attention is necessary to the design of Characteristic Space because the result becomes the one that differs by the different design. In this case, it is sufficient to use the raw data itself as Characteristic Space because the raw LiDAR data is the position information. It is often used to search the neighbor of a point applying the Cluster Analysis. How it defines the neighbor is the key to the Cluster Analysis, Voronoi Neighbor of the XY plane is adopted here. There is a possibility that it becomes the neighbor in XYZ space even the points have a height difference because the raw LiDAR data is random. Also, there is the reason that the Voronoi segmentation is used as the building shape mentioned in Section 2.6. Furthermore, it needs to be classified uniquely without the effect of the input order of the point data. In other words,

The comparative function is derived from a height difference in Z value or an inclination form. The threshold value T is defined as a function. For example,

As taking another neighbor, a point within a radius on XY plane is conceivable. Therefore, Cluster Analysis of raw LiDAR data is carried out with the following rule.

2.2 Surface Determination Part 1

The method that gives the recognition of the terrain and/or non-surface objects to the raw data that has been clustered by the height characteristic mentioned in Section 2.1 is discussed here. Each cluster is determined whether it is terrain or non-surface objects depending on the land model. First of all, it is known by intuition that the terrain has many number of the clusters. Therefore, the cluster that has over D~S1 [points] is determined as the terrain when the average density of the raw data is D [points/m2] and area of the biggest building is S1 [m2]. Next, elimination of cars and roadside trees from the raw data is considered. These areas are smaller than the area of the smallest building. Therefore the cluster that has below D~S2 [points] determines cars or roadside trees when the area of the smallest building is S2 [m2] (Table 1).

2.3 Surface Determination Part 2

The determination method of the terrain and/or non-surface objects to group B that has been classified in Section 2 is discussed here. Group B is determined whether it is the terrain or not depending on the basis of group A that has been already determined as terrain. There are various methods in Cluster Analysis and the problem is "the distance function" and "how to take distance between clusters". In this case, "Euclid Distance of the XY plane" and "Nearest Neighbor Method" are adopted for the distance function and methodology.

2.4 Return Information Analysis

There is a problem that the last pulse is not necessarily reaching the surface, although the return information is effective in forest areas as mentioned in the beginning of Section 2. Therefore, the determination whether the last pulse is the terrain or not in forest areas is necessary. It has been explained in Section 2.3 about the determination whether objects are the terrain or not. Thereupon, the last pulse is also determined by adding itself to group B.

2.5 Generating DTM

Generally the DTM is generated from interpolation with group A that is judged by Section 2.2 and group B1 judged by Section 2.3 as terrain. It is still a problem, however, to express the terrain where the building exists. There is a possibility that such area becomes sloping ground by interpolation method. To avoid it, the terrain of this area is estimated from the neighbor terrain of group B2. According to this method, as another problems, it may have been determined as the terrain even the bridge and overpass etc. In any case, these problems are influenced by the definition of the DTM.

2.6 Generating Building Shapes

A building shape is generated from group B2 that is judged as non-surface objects in Section 2.3. A boundary line is calculated by connecting each of Voronoi region by group B2. However there are several problems. If there is no data available between the buildings, it is recognized as one building. On the other hand, it becomes a separated building if there is some height difference or rooftop structure even it is one building. In addition, a building height is calculated as follows, and how to determine the terrain height is the problem as expressed in section 2.5.

These problems are also influenced by the definition of the building shape and its height.

2.7 Tree Extraction

It is assumed that the point other than the terrain is a tree in the forest areas. The calculation of individual tree height is difficult because it is indistinct where the laser pulse hits and be reflecting on the tree. However it is possible to calculate the volume of trees from the height and Voronoi region and the point density.

3. Test and Result

The test was carried out by applying the algorithm as mentioned in Section 2.

3.1 Urban Areas

The data shown in Table 2 was used for extracting DSM in Urban Areas. Figure 1 shows the data processing flow. First of all, inspection of the terrain is done. Figure 5 shows the sectional combination of the raw data (Figure 2), the interpolated raw data (Figure 3) and the DTM (Figure 4). The points that are catching the terrain generate the DTM. Also, the abnormal values could be removed as another effect. The abnormal value occurs with diffused reflection of water area, window, etc.

Furthermore, the buildings that are generated with this processing (Figure 7) were compared with the buildings by map digitizing (Figure 6). The several outlines of the buildings match in suburb areas because the interval of building is wide. But the result was inadequate in urban areas (Figure 7). However, the fact that the point could be generated automatically and building could be extracted without the influence of cars and roadside trees leads to the next step.

Finally, the DSM (Figure 8) was generated from the DTM (Figure 4) and building shapes (Figure 7). Comparing this to the Figure 3, it is evident that terrain and buildings are separated completely and cars and roadside trees are eliminated.

3.2 Forest Areas

The data shown in Table 3 was utilized to extract DSM in Forest Areas. Figure 9 shows the data processing flow. At first, the inspection of the terrain is conducted. Figure 13 shows the sectional combination of the raw data (Figure 10), the interpolated raw data (Figure 11) and the DTM (Figure 12). The terrain could be extracted even in the area of dense trees by using last pulse. And the last pulse that is not catching the terrain could be eliminated.

Figure 14 shows the combination of the contour generated by this processing and by stereo plotting. Although it matches well in non-dense area, the contour line does not match in dense area. This is the influence of the difference of the processing method in the dense trees area. The digitized contour by operator was drawn with prediction of surface by visual inspection. As the acquisition of a last pulse becomes more difficult in this processing, the influences also become bigger. In any case, there is a limit to complete terrain extraction in the dense tree area. So it should be supplemented by other method.

Further, trees are extracted by deducting the terrain from the raw data (Figure 15). This data is not the exact tree canopy because it is indistinct where the laser pulse hits and remains reflecting on the trees. However, the fact that the LiDAR point data could be generated automatically and trees could be extracted certainly leads to the next steps.

Finally, the DSM (Figure 16) was generated from the DTM (Figure 12) and trees (Figure 15). Comparing this to the Figure 11, it is also evident that terrain and trees are separated.

4.Conclusions

By this proposed method, automatic computer DTM extraction from LiDAR data was completed as a manual classification with visual inspection. The extraction degree of this method is higher than the edge extraction method of traditional image processing. However, the satisfying result could not be obtained in the field of the shape extraction of non-surface objects. The reason is to face the difficulty for the extraction of shape from buildings that have various rooftop conditions in urban areas and also the indistinctness in forest areas where the laser pulse hits and remains reflecting on the tree. The extraction of non-surface objects from the LiDAR point data is difficult even if the judgment will be done visually. This fact means that there is a limit in the field of DSM extraction from LiDAR data that is composed of the point data.

Some problems, e.g., automatic recognition of urban areas and forest areas even in the terrain extraction still exist. The urban area where building exists is fully developed and the model differs from the forest. It is necessary to develop method to divide urban areas and forest areas in order to raise the accuracy of the surface extraction. To increase the accuracy of surface extraction, a proper processing method should be adopted depending on the areal characteristics.

The image information is indispensable to solve these problems. High resolution multiple spectrum airborne sensor has been developed in recent years. The next target of research is the development of high accuracy DSM extraction depending on the combination of the multiple spectrum image by the CCD line sensor and LiDAR data.

Acknowledgment

The author would like to thank all for their helpful advice on the writing of this paper.

References

1) Raul Campos-Marquetti, Robert Kletzli, John Nipper, Scott Paulsen and Steve Scharf, (2000). Digital Terrain Mapping of Bernalillo County, New Mexico using Digital Orthophotography and Airborne LIDAR. URL: http://gis.Esri.com/library/userconf/proc00/professional/papers/PAP726/p726.html.

2) Kinji Mizuno, (1996). Analysis of Multivariate data, Asakurashoten, Tokyo.

3) Takeshi Agui and Tomoharu Nagao, (2000). Introduction to Image Processing using Programming Language C, Shokisya, Tokyo.

4) Shinya Ishikawa, (2001). URL: http://www5.ocn.ne.jp/~shinya91/.


Masaomi Okagawa
Geomatics Team, GIS Institute
Pasco Corporation
E-mail: masaomi_okagawa@pasco.co.jp