Rule-based Learning Techniques to Derive Automated Digital Terrain Model Using Airborne LiDAR Data

DOI: https://doi.org/10.52939/ijg.v18i6.2459

1. Introduction

The surface shape, shown by the Digital Terrain Model (DTM), displays the earth’s surface using points with specific coordinates [1] in the raster format, which is one of the main products of photogrammetry and remote sensing [2] and [3] and the basis for many applied projects [4] [5] and [6]. It has always attracted the attention of many experts because the accuracy of the DTM can dramatically change the research results [7]. However, many research works that use DTM as an information layer, employ commercial software [8], the majority of which use the data without applying necessary filters, although filtering contributes substantially to creating an accurate DTM [9] and [10]. In addition, filtering LiDAR data needs time and specialty, which depends on the experts' experience [11] and [12] and is not always possible for many executive and research projects. The importance of this issue becomes more apparent when carefully reviewing the research by Klapste et al., [13] on DTMs generated by some software. Although the area studied by them did not have man-made structures, and the experts had previously visited the points, the RMSE results showed values in the range of 0.13 to 0.19m. Therefore it is essential to strictly investigate the accuracy of the generated DTMs by the common software in the urban areas, especially where no initial review of the LiDAR points has taken place.

Furthermore, the need for an accurate DTM in studies based on extensive details such as surface water and flood management [14] is an undeniable fact. However, from the theoretical point of view, LiDAR can create a digital elevation model of the earth’s surface with high speed and accuracy compared to the traditional methods [15] [16] and [17]. Yet, the extraction of the digital elevation model from urban areas and complex and uneven environments is challenging, even using modern methods.

Hence, the accurate DTM is only obtained when the non-ground points are properly filtered [18] and [19] and correctly interpolated. Despite all studies conducted, building an accurate DTM has always been a serious challenge in applied research [20] and [21]. The methods for accurate generation of DTM are applicable only when the experts can use them like the available commercial software without the need for specific knowledge.

1.1 Literature Review

Perhaps filtering based on a specific window can be used as one of the oldest DTM development methods that have evolved. Kilian et al., [22] and Zhang et al., [23] used this method, indicating that if the correct information on the region is provided, the resulting DTM has adequate accuracy. The slope-based method has been also used to separate the ground and non-ground points. This method initially considers a point as a ground point, which serves as the basis for the measurement of other points [24] and [25].

Sohn and Dowman [26] proposed a method, which first created a network with several initial points, and then a decision was made based on the angle of each point concerning the surrounding points to omit or keep the point as a non-ground or ground point, respectively. To improve this method, the initial points for triangulation were selected from the corners of the region, and the remaining points were classified into the upper and lower groups of initial triangulation. Afterward, the conditions of triangulation angles were studied separately for each group.

In 2009, Meng et al., [18] attempted to compare the properties of each point with the four points around it and created a filter based on four different types. Although this filter had satisfactory results, interpolation based on the independent value of each point, especially in large buildings, was not satisfying due to the drastic height variations on the building boundary.

Baligh Jahromi et al., [27] used an artificial neural network to omit the non-ground points. This method is based on the height difference between the first and last reflections. However, it cannot be used to develop DTM with fast processing due to the multiplicity of the constraints on the angles defined and the need for sample data. On the other hand, steep slopes are a limiting factor in this method, leading to numerous challenges in the practical application of this method in operational programs.

In 2016, Zhang et al., [28] conducted a study of DTM extraction from DSM and stated that it is substantially important to separate non-ground points from the ground points for the automatic extraction of DTM from LIDAR data. This study sought to extract DTM from DSM with low resolution, using a two-step filter. The experiments on the resulting DTM showed significant effectiveness of the method for data with spatial resolution of 10 to 30 meters. However, the algorithm can be only used for large areas to provide an overall view and cannot be used for urban areas that call for high accuracy.

Wojciech [32] carried out a study on building DTM using the kriging method and tried to examine the usual interpolation methods to create an accurate DTM. They examined the inverse distance weighing, the nearest neighbor, the moving average, and kriging algorithms. This study aimed to analyze the characteristics of the results of each of the existing interpolation methods to considerably reduce the computation time while maintaining the highest accuracy level in the created model. The experimental results highlighted significant advantages of the kriging technique over the other methods [29] [30] and [31].

1.2 Objective and Research Gap

Generating accurate and applicable DTMs remains a serious but essential challenge [10] [31] [33] and [34]. Although all the previous studies have sought to generate accurate DTMs, due to the complexity of the used filters [11] and [32] in a section of research process, they have taken advantage of one stage of the specialist supervision. This would result in two problems: 1) when the DTM layer generation is not the primary objective, but researchers utilize DTM as one of the necessary layers of information, it is not possible to implement essential modifications through field visits or taking advantage of experienced experts. Thus, they prefer to use raw data and commonly available software. 2) As the factors, which the experts should examine for the generation of DTM, are numerous and related to the human experience, certainly there will be errors during the investigation process. Therefore, the methods presented in the research works are only operative when the proposed algorithm can generate the required DTMs automatically without the need for further data or expert information. In reality, the most authentic method, which is the manual method [17], requires extensive time and costs, making its application impossible in the research and executive projects.

Therefore, an important point to note is that the DTM has applications in a broad spectrum of research works and operational projects. In many projects, DTM is just one of the layers used, which can significantly affect the results.

Previous research works have largely ignored the fact that although a complex, time-intensive method can lead to an accurate DTM, this process is not practically feasible for other research works and operational projects as it demands for long processing times and an expert who is familiar with the generation process of DTM. Accordingly, a clear research gap in the existing literature is the failure to pay adequate attention to the practical applicability of the proposed DTM production methods to ensure their usability for all researchers and experts who have no expertise in DTM generation.

Therefore, creating a fully automatic structure is one of the main goals of this research. For this purpose, the LiDAR data are implemented because LiDAR data provide the possibility to transfer geographical location and elevation characteristics and much other information about the land objects in the form of point clouds. On the other hand, it is possible to convert LiDAR data into a multidimensional matrix, which subsequently helps with using the programming languages in the analysis of the above mentioned data. This research has the Python programming language to write a program which initially accepts the LiDAR data file, executes a number of rule-based algorithms, and finally stores the DTM file with the raster format in the exit path.

2. Materials and Methods

2.1 Study Area

The study region with a surface area of 3000ha is located in the west of Malaysia, southeast of Selangor Province. The coordinates of the aforementioned region in the UTM coordinates system in the region N47 are as follows: minimum x = 800000, maximum X = 805000, minimum Y = 328000, and maximum Y = 334000 (Figure 1). The study region is divided into three major parts concerning texture and land use: 1) Areas including one-story buildings and lands covered with grass and scattered trees, highly suitable for research in rural areas; 2) Dense urban areas with high-rise buildings, urban drains, power transmission lines, and low vegetation, suitable for research in urban areas; and 3) Areas without urban texture, covered with grass, trees, canals, and scattered lakes, and suitable for research in non-urban regions. Since the study area was extensive, and the suggested algorithm's findings required repeated field visits, access to the study locations was critical. As a result, the research region was separated into three sections: urban, rural, and without structures, each sector including accessible samples of 2 hectares.

Finally, the suggested algorithm assessed 301 hectares of the region, of which 136 hectares were urban areas, 88 hectares belonged to rural areas, and 77 hectares were without buildings.

Figure 1: Location of study area

This research used LiDAR data for all processes and RGB images for visual inspection of manual DTM construction and preliminary checking of results. LiDAR point cloud data and orthorectified, used in this research, were acquired in 2015 over Universiti Putra Malaysia and part of the Serdang area by Ground Data Solution Bhd. The raw LiDAR dataset was collected using the Riegel scanner on-board EC-120 Helicopter hovered at an average height above sea level of 600m above the terrain surface. The orthorectified (RGB color image) was captured with the Canon EOS5D Mark III camera with a focal length of 35mm, a horizontal and vertical resolution of 72Dpi, and an exposure time of 1/2500sec. The density of LiDAR data was about 8 points per square meter. This research performed data processing in five steps, all of which were taken automatically in an algorithm written using the Python language.

2.2 Separation of the Last Reflection

Each point in the LiDAR data carries several descriptive data. In the data used in this study, the initial separation is based on the first and last reflections, and the points in the LiDAR data carry a factor known as the ‘Classification factor’. Therefore, x, y, and z determine the position of the point in space, and the class factor can determine the first or last reflection. At this stage, the point clouds are converted into a Numpy array, the points with the last reflection are identified, and the other points are removed using the classification factor in the array.

2.3 Calculating the Effective Distance

The designed algorithm compares each point with the surrounding points at a specific distance. More importantly, the effective distance is different for each data, and the above algorithm is designed to automatically calculate the effective distance in each row and between the rows for each different data. The effective distance in each row is larger than the average distance between two points and smaller than twice the average distance between two points, always including at least one point. The effective distance between the rows refers to a distance larger than the average distance between two rows and smaller than twice the average distance between two rows. The average distance between the points in each row is calculated based on Equation 1, and the effective distance in each row is calculated based on Equation 2. However, it is worth noting that this distance is calculated based on the horizontal distance without deviation from the main directions. To operationalize the effective distance, the value that has to be added to the x value of each point is calculated based on Equation 3. D2 is the value whose elimination or addition to the value X of each point will create a range that contains at least two other points except the target point.

Equation 1

Equation 2

Equation 3

Equation 4

Equation 5

Equation 6

Equation 7

Where:

x	Value of x axis
y	Value of y axis
i	Number of points
n	Number of points available in LiDAR data
j	Number of rows
m	Number of rows available in LiDAR data
d1	Average distance between the points in row
d2	Effective distance in each row in algorithm
d3	Average distance between rows algorithm
d’	Effective distance between points
d’’	Effective distance between rows

The effective distance for each row is also obtained based on the average distance between the rows. Therefore, the line equation for the row is obtained according to Equation 4, and the distance from each point to the line is calculated based on Equation 5, followed by the calculation of the effective distance for that point using Equation 6. The resulting distance actually determines the vertical distance from each point. To operationalize the results and use them in the designed algorithm, the resulting effective distance has to be converted into a number to be added to and subtracted from the value of y of each point. Therefore, the resulting value is converted to the distance value on the y axis using Equation 7. Based on the calculations above, a range is identified for each point based on the minimum and maximum values calculated on the x and y axes. As a result, at least 4 points are placed around the target point.

2.4 Error Filtering

Errors in LiDAR data are classified into two categories: Type I includes ground points regarded as non-ground points, and type II includes non-ground points accepted as ground points [35]. In urban areas, numerous high-density man-made objects can cause type II errors. Therefore, the filtering stage is substantially important in LiDAR data and can be performed both manually and automatically [4]. The current research has used the nearest neighbor filtering based on the effective distance because this filter can decide according to the features of the data existing in the neighborhood of the target point, and this is fully in line with the selection of a point for the grounds and non-grounds. Furthermore, the research shows that this filter is effective at noises removal [36].

2.4.1 Nearest Neighbor Filter (K-NN)

Here, the effective distance is the range required to analyze the data features. The logic behind this idea is that if the selected point is part of the ground points, most of the points around it fall in the same category, and if a point is a noise, most of its neighboring points do not fall in this category. Therefore, the K-NN filter can remove many noises and the ‘roof error’ is one of the errors removed with this filter.

2.4.2 Roof error

Laser can enter some objects due to their characteristics and lead to the initial and last reflections. This naturally applies to trees, allowing for the recognition of vegetation. However, in urban areas, this phenomenon occurs in some ceilings through which laser can pass and create what is known as the ‘roof error’. These types of structures are not exclusively used in greenhouses with transparent ceilings because urban grounds have various structures, and this error is widely present in a part of the existing structures. Considering the logic used in the neighboring filter based on the effective distance, the K-NN filter can remove the roof error. Pseudo-code as illustrated in Figure 2 is a part of codes written for the K-NN filter automation in the DTM building stage.

Figure 2: Pseudo-code for K-NN filter automation

2.5 Creating a Modified Network

After filtering the nearest neighbor, many locations lacking data, known as the ‘empty spaces’, are created. These empty spaces are classified into two categories. The first category includes small spaces created due to the removal of one or two points, and the second category includes large spaces created to remove the structures. In this research, ‘large spaces’ refer to spaces in which there is no remaining point within the effective window range. A new grid of points has been used to fill the empty spaces. The distance between points in the new grid is assumed equal to the minimum distance between the points in LiDAR data and rated based on the elevation value of the points in LiDAR data. Hence, the weighted averages of all the existing points in the LiDAR data, located within the effective window with the center of target point in the new grid of points, are calculated. Therefore, the value of the target point is calculated based on the distance of each LiDAR point, which is in effective window from the target point according to Equation 8. This process is repeated for all the points in the new grid. The value of existing points in the large empty spaces is filled with -999.

Equation 8

where:

Z(So)	Elevation at target point
Z(Si)	elevation at point i
i	considered point number at effective distance
n	number of points in effective window
k	impact coefficient of target point distance, which is considered here as 2
so	target point
si	point i in effective distance

2.6 Interpolation

This stage targets the large empty spaces filled with value -999. First, the boundary of every object is detected because the margins of structures are significantly important in urban planning. This is because a structure is built as a vertical object into the ground, and only the structure boundary is in line with the general earth slope, which is considered in environment planning such as the design of the urban sewage network. Hence, if a structure is removed, the space occupied by the structure has to be filled as a smooth surface in line with the general earth slope.

In other words, if the boundary points maintain their effects in the object interpolation, two types of error occur, including stair and surface interpolation errors.

2.6.1 Stair error

Stairs are objects sometimes implemented to connect land to buildings. They are usually seen only in the input section and are not present at all the borders of the building (Figure 3(a)). However, in the design of the urban drains, the general slope of the ground is often followed, which is the same slope of the building border (Figure 3(b)). Hence, the space under the stair is often used to continue the drains (Figure 3(c)). In other words, the stairs connected to the buildings do not act as a barrier against running water. Therefore, when filling the empty spaces, if interpolation is done normally and considering all the existing points in the area, the impact of steps connected to buildings is observed within the buildings limits. In this case, the spaces within the buildings would turn into wavy surfaces instead of being flat (which follow the elevations at the borders of the building). This type of surface could not be a good DTM for the urban area. Here, this type of error, produced due to any object that connects land to a part of the building, is called the “Stair error”.

2.6.2 Surface interpolation error

All the existing interpolation algorithms are designed based on data continuity, and the continuous data maintains its effect in the new data up to a certain distance. From there, water flow follows the overall slope of the structure. Therefore, an algorithm has to be used to create a smooth part after interpolation. As a result, the designed algorithm first modifies the boundary of the large empty spaces based on the average height of each side, and interpolation is performed correctly within the range. In this case, the effect of stairs or small structures, which play the role of connectors, is eliminated, and the impact of the height of the sides is not stretched along the structure.

2.7 Create DTM

At this stage, the new point’s network, edited in the previous steps based on the coordinate system of LiDAR, is converted to DTM. Therefore, the DTM is created in a fully automatic process from LiDAR data and saved in the TIFF format in the output path.

2.8 Accuracy Assessment

Since the manual DTM has good accuracy [34], the DTM of the area was first manually generated in this section, using the LiDAR point clouds, and both DTMs were compared based on statistical calculations and field visits.

2.8.1 Statistical calculations

A database was created, where each point had two values of DTMs, manual and automatic. Descriptive statistics were used to evaluate the height and slope of the two DTMs. The minimum, maximum, and average height and slope factors of the two DTMs were investigated, and their standard deviation was calculated based on Equation 9. In the next step, the difference between the two layers was examined via the t-test. Finally, the two DTMs were compared to each other based on RMSE because it is an appropriate criterion to examine the amount of error [33] [37] [38], and [39].

Equation 9

Where:

N	Number of available database points
I	Variable number under review
xi	i-th variable under review
µ	The mean of the data under review

2.8.2 Field visit

Considering that some differences in the output layers required a field study, a change map for the two DTMs relative to each other was generated, and the coordinates of the changes were extracted from the change map. Next, using the RGB image, the object that caused two different DTMs was determined, and the automatic and manual DTMs were compared in the ArcScene environment in a three-dimensional form. In the next step, using GPS, the mentioned object was found in the study area, and the difference between the two DTMs were subject to field visits.

3. Results

3.1 The First Study Area

The first study area has one-story buildings, narrow roads, and relatively dense tree vegetation that are sometimes very close to the buildings or cover them. Furthermore, the construction of buildings does not follow a particular order, and temporary structures are abundantly scattered everywhere. This region has all the characteristics of rural areas. In Figure 4(a), the building shown with the letter D is considered an object using the algorithm designed in this research, the boundary of the given object is modified based on the designed algorithm, and then interpolation is performed. As a result, the canal on the left side of the building has remained a real canal and does not affect the building environment. However, in Figure 4(b), which is associated with the manually built DTM, the height of the existing canal directly affects the location of the building due to the lack of the mentioned mechanism. As a result, this canal appears as a series of continuous pits with separating hills, while the building area looks like a pit. Figure 4 clearly shows that the removal of the ‘stair error’ and the ‘surface interpolation error’ is necessary to correctly display the building locations and canals. A comparison of the elevation and slope characteristics of the corresponding points in both DTMs shows no significant differences between the elevation factors (Table 1). However, there were evident differences in the slope factor. Next, the value of each DTM was subjected to the independent samples t-test (Table 2), revealing a significant difference between the two DTMs. Although the above test indicated a difference between the two DTMs, the t-test could not show the value of this difference. This issue needs further investigation, considering the correlation of 95% between the two DTMs, the standard deviation value of 0.001m, and the error range of 0.013-0.018m, as it might be impossible to use the difference for specific scales.

Hence, the RMSE value of the DTM layer created automatically was calculated considering the DTM calculated manually. According to the calculations, the RMSE was .05m in the rural areas, indicating that the average difference between automatically and manually calculated DTMs was 5cm. The used maps in the urban and rural planning had the scale of 1:2000 [40]; therefore, the error of 5cm is totally acceptable.

3.2 Second Study Area

This area is a prominent example of an urban area. The presence of high-rise buildings, numerous cars, wide streets, and concrete structures such as curbs and short trees compared to the buildings are the characteristics of the second study area, which was examined as an urban area in this study, with the selected samples covering an area of 136 hectares. Field visits showed that many buildings in the second study area have a longitudinal slope and are built as stepped buildings on a slope. However, the manually built DTM is entirely uniform and does not consider the changes made to the earth’s surface. Figure 5 shows the location of the investigated samples within the study area while also demonstrating a residential house to better illustrate the longitudinal and lateral slopes and highlighting the difference between the two DTMs.

Figure 3: A view of the design procedure of urban drains with respect to buildings

Figure 4: Constructed DTMs. (a) automatically generated DEM (b)Manually generated DEM

Table 1: Slope and elevation specifications of constructed DTMs in first study area

Table 2: Comparison of two DTMs based on independent-sample t-test in first study area

Figure 5: Longitudinal and cross section of manually and automatically constructed DTMs

As clearly seen in Figure 5, the ground at the bottom of the building is slightly hollow in the longitudinal cross-section of the manually built DTM, and it definitely cannot be a suitable model for urban planning. However, the DTM built based on the proposed algorithm fully restores the changed ground shape and shows the increased height at the building downstream. This is manifested in another form in the transverse cross-section of the building.

Here, the comparison of the elevation factors of the two DTMs indicates a slight difference between the study factors. However, an analysis of the existing slope difference revealed a significant difference (Table 3). Seemingly, as the number of urban structures increases, the difference in the average slope increases, and thus the sudden slope changes in urban regions are considerably higher than the rural regions. The statistical test also indicates the difference between two DTMs (Table 4); however, the value of difference is not clear. The correlation of 95% between the two DTMs, standard deviation of 0.007, and an error range of 0.047-.073 confirm the slight difference between the two DTMs. Hence, the RMSE of the DTM layer constructed automatically was also calculated considering the manually-constructed DTM. The calculations show the RMSE value of 0.2m in the urban areas, highlighting that the average difference between the DTMs constructed automatically and manually is 20cm, which is totally acceptable considering the required scale in the urban and rural planning.

3.3 Third Study Area

This area lacks man-made buildings and structures and is covered with grass and scattered trees. Hence, this region is highly suitable to analyze the mechanisms used to build the DTM manually and automatically. The results indicate that as the number of urban structures grows, the difference between the manual and automatic methods is further highlighted. Therefore, if the difference between the two methods used is based on man-made structures, the first hypothesis is that in areas lacking man-made structures, there should be no differences between the two DTMs built.

Table 3: Slope and elevation specifications of constructed DTMs in second study area

Table 4: Comparison of two DTMs based on independent-sample t-test in second study area

Table 5: Slope and elevation specifications of constructed DTMs in third study area

Table 6: Comparison of two DTMs based on independent-sample t-test in third study area

Figure 6(a) shows the DTM built using the proposed algorithm, and Figure 6(b) shows the DTM built manually. The results of the field studies confirm the primary hypothesis, according to which there is not a significant difference between the two DTMs. A comparison of the elevation characteristics of the two DTMs also shows many similarities, as confirmed through the comparison of the slope characteristics (Table 5). This is also confirmed by the t-test, and the numbers in the two DTMs linked to one coordinate are not significantly different (Table 6). On the other hand, the RMSE=0.04m exhibits an error equal to 4cm, which is a totally acceptable precision.

Figure 6: DTMs constructed in third study area. (a) automatically generated DEM

(b) Manually generated DEM