# The Estimation of Tree Height Based on LiDAR Data and QuickBird Imagery

Wei Su, Rui Liu, Ting Liu, Jianxi Huang, Xiaodong Zhang, Junming Liu

► To cite this version:

Wei Su, Rui Liu, Ting Liu, Jianxi Huang, Xiaodong Zhang, et al.. The Estimation of Tree Height Based on LiDAR Data and QuickBird Imagery. 6th Computer and Computing Technologies in Agriculture (CCTA), Oct 2012, Zhangjiajie, China. pp.472-482, 10.1007/978-3-642-36137-1_55. hal-01348265

HAL Id: hal-01348265

https://inria.hal.science/hal-01348265

Submitted on 22 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# The estimation of tree height based on LiDAR data and QuickBird imagery

Wei Su¹,², Rui Liu², Ting Liu², Jianxi Huang², Xiaodong Zhang², Junming Liu²

¹Key Laboratory of Agricultural Information Technology, Ministry of Agriculture, Beijing, 100081, China;

²College of Information and Electrical Engineering, China Agricultural University, Beijing 100083

**Abstract:** The estimation of tree height is advanced following the development of LiDAR technique. The estimation model of tree height considering suppressed trees is developed in order to extract tree height accurately using LiDAR data. Filtered LiDAR data and Quickbird imagery are segmented using watershed segmentation method based on mathematical morphology to get the boundary of trees. And the highest point in each canopy object is used to estimate tree height. Weibull distribution is used to estimate height distribution of the suppressed trees. The experiment results indicate that the watershed segmentation method based on mathematical morphology is an effective method to extract the boundary of trees. And the R² between the tree height estimated using estimation model of tree height considering suppressed trees and the tree height measured by field work is 0.93.

**Keywords:** estimation of tree height, LiDAR, watershed segmentation, Weibull function

## 1 Introduction

LiDAR (Light Detection And Ranging) is a remote sensing technology that can be used to extract detailed 3D structure information. LiDAR is an active system which emits high frequency bursts of laser pulses. So they are able to penetrate relatively dense vegetation canopies, and can be used to get the 3D structure information of trees even for suppressed trees.

There are many researchers explore how to estimate tree height using LiDAR data in recent years, most of them do this using CHM (Canopy Height Model). This CHM is the difference between DSM (Digital Surface Model) and DTM (Digital Terrain Model), which is the height of non-ground objects such as trees, buildings, vehicles etc. Liu et al. develop smoothed CHM to estimate tree height. Firstly, they smooth CHM using Gaussian smoothing method. Then, they develop a double tangent crown edge recognition (DTCER) algorithm to detect the canopy boundary, and the local maximum search algorithm with fixed or variable window was used to recognize crown tops. Persson et al. create DCM (Digital Canopy Model) from subtracting DTM from DSM, then which is smoothed with different scales. A parabolic surface was fitted to the elevation data to determine which scale to choose for different parts of the image, and the height and crown diameter are estimated for the identified trees [7]. Popescu et al. use regression models and cross-validation to estimate tree height based on LiDAR data. They found that estimating tree height for deciduous plots gave superior results without calibrating the search window size based on forest type [8]. Zimble et al. search potential canopy top points using a moving radius circular search window with 0.9m radius, and adjoining pixels that are higher than 85% of their
neighbors are combined into clumps that are assumed to include the peaks of tree crowns[9]. A common ground of all of these researches is that there should be some prior knowledge about their study area which is the basis of determining initial parameters used in model. Therefore, there are some difficulties in extending those methods to other areas.

We develop an estimation model of tree height considering suppressed trees in this study. The watershed segmentation algorithm is used to segment tree canopy objects, the highest point in each tree canopy object is searched using a moving window size, and the height of suppressed trees is estimated using Weibull function. The rest of the paper is organised as follows: Section 2 describes the study area and data source. Section 3 presents the methodology of the estimation model of tree height considering suppressed trees. And the results, discussion and conclusions are presented in Section 4 and Section 5.

## 2 Study area and data sources

The study area is an urban area (Kuala Lumpur City Centre, from 101°40'46" to 101°43'25"E, 3°8'00" to 3°10'39"N) located midway along the west coast of Peninsular Malaysia, at the confluence of the Klang and Gombak rivers. This region is characterized by rolling topography and flourishing vegetation, and the height of high trees range from 20m to 30m. A subset area covered with dense high trees in Kuala Lumpur City Centre is selected as study area in this study. And the QuickBird imagery and CHM image produced using LiDAR data of this subset are as Figure 1 (a) and (b) respectively, which are resample to 1m spatial resolution.

<figure><img src="image_2.png" /><figcaption>Fig.1. Study area (subset of Kuala Lumpur City Centre, Malaysia). (a) Quickbird image, (b) CHM image</figcaption></figure>

## 3 Methodology

The developed estimation model of tree height considering suppressed trees works on CHM and Quickbird image. The ground objects such as trees, buildings, roads, are segmented based on CHM image using watershed segmentation based on mathematical morphology method. The spectral features coupling with height information of image objects are used to identify canopy objects. So the boundary of
canopy are acquired. Searching the highest points in each object, we can estimate the height of trees. At last, we use Weibull function to simulate and estimate the height of suppressed trees. And there are 4 components in this algorithm: CHM image smoothing, image segmentation using watershed segmentation method based on mathematical morphology, extracting the tree vertex points and estimating the height of canopy, estimating the height of suppressed trees using Weibull function. This flow chart of the whole algorithm is as follows:

<figure><img src="image_3.png" /><figcaption>Fig. 2. Flowchart of estimation model of tree height considering suppressed trees</figcaption></figure>

## 3.1 Image smoothing

Subtracting DEM from DSM, we can get CHM of canopy. However, for the irregularity of LiDAR data and the twigs of the canopy surface, there will be false tree vertex if we use unsmoothed CHM image directly to identify the tree vertex and crown boundary. So we smooth the CHM image in order to remove these noise before identifying tree vertex and crown boundary[10]. However, this smoothing will weaken the difference between tree vertex and around pixels, which will decrease the ability to identify the tree vertices. So he Gaussian filter method is used to smooth CHM image. The filtering kernel of Gaussian filter is as follows:

$$g(x, y, \sigma) = \frac{1}{2\pi\sigma^2} e^{-(\frac{x^2+y^2}{2\sigma^2})} \quad (1)$$
Where, x and y are the distant to the center of smoothing nuclear, σ is used to adjust the parameters of the Gaussian function. Figure 2 shows a schematic diagram of a 3×3 filtering kernel of Gaussian smoothing algorithm.

<figure><img src="image_4.png" /><figcaption>Fig. 3. 3×3 filtering kernel of Gaussian smoothing algorithm</figcaption></figure>

Moving this filtering kernel, we smooth CHM image. For the pixel p(x,y) located in the centre of filtering kernel, its DN value is computed as:

$$p(X,Y) = \sum_{i=-d/2+0.5}^{d/2-0.5} \sum_{j=-d/2+0.5}^{d/2-0.5} p(X+i,Y+j) \times g(x+i, y+j, \sigma) \quad (2)$$

X and Y are coordinates of the pixel in the image, d is the size of the filtering template which is an odd number such as 3, 5, 7 etc., i and j are the distant to the studied pixel.

## 3.2 Watershed segmentation based on mathematical morphology

### 3.2.1 Principle of watershed segmentation

The watershed segmentation algorithm is a robust image segmentation tool, which can be used to extract interested objects in the image effectively. Vincent's simulate immersion algorithm is used in this study [12]. The pixels in gradient image are classified as four types during the watershed segmentation: marked pixel, waterline pixel, unlabelled pixel and classified pixel. The watershed is finished when all pixels are classified pixels or waterline pixels. There are two steps in Vincent's method: sorting and immersion.

1. Sorting: ascending the pixel value in line with the DN value of pixel, we can get a sorted pixel matrix. In the process of gradually submerged, we can call the pixels which need to be processed. But we needn't process all pixels in every loop, which can speed up the calculation.

2. Immersion: this immersion is performed in ascending order to call each pixel based on sorted gradient image. Starting from the minimum DN value in the whole gradient image, immersion is done one pixel by one pixel. This loop queue is used to extend labeled catchment basins. Based on the distance to the neighborhood, we label finished watershed.
<figure><img src="image_5.png" /><figcaption>Fig.4. Flowchart of Vincent's segmentation algorithm</figcaption></figure>

The study area is a tropical rain forest area, where there are dense and high tropical trees. This results in a difficulty in canopy segmentation because there is overlapping between one canopy and another. Watershed segmentation method is an effective way to extract canopy objects because which can take advantage of the height difference between trees fully. In addition, the watershed segmentation method is a method without any parameters, which is suitable in many situations. So we can get an ideal image segmentation result without prior knowledge. However, there will be a minor gradient change if there are two trees which is located near and the canopy top is relatively flat. In this case, it is difficult to identify these tree objects using watershed segmentation based on gradient. Fortunately, watershed based on the distribution of height is used to do further segmentation. The watershed segmentation algorithm based on the distribution of height is as follows: searching the maximum points in each objects produced from previous segmentation using fixed window, the further segmentation will be done if there are two or more maximum points in one objects.

Waterline pixel in the watershed segmentation is the intersection of two or more tree objects. The classification of these points is done by the inclination angle between this point and crown vertex which the studied point is possible belonged to. And the studied point should be classified as the tree objects which have a larger inclination angle. For an example in Figure 4, there is an intersection between tree A and tree B. r1 and r2 are the inclination angles between this intersection and tree A and tree B. And this point should be classified as tree A, otherwise, classified as tree B.

<figure><img src="image_6.png" /><figcaption>Fig.5. Schematic diagram of ownership of the intersection point</figcaption></figure>
### 3.2.2 Merging of over-segmentation

There is over-segmentation in watershed segmentation sometimes, especially when there are larger branches highlighted in studied objects. So we should merge these over-segmented image objects. The Full Lambda Schedule algorithm of Robinson et al. is used in this study[13], which is as follows:

$$t_{i, j} = \frac{|o_i| \times |o_j|}{|o_i| + |o_j|} \times \|u_i - u_j\|^2 / \text{length}(\alpha(o_i, o_j)) \quad (8)$$

$|o_i|$ and $|o_j|$ are the area of image object i and j respectively, $u_i$ and $u_j$ are the mean DN value of image objects i and j, $\|u_i - u_j\|^2$ is the Euclidean distance between the image object i and the image object j, $\text{length}(\alpha(o_i, o_j))$ is the length of the public side of the image object i and the image object j. $t_{i, j}$ is computed using Formula (8), we will merge image object i and image object j when $t_{i, j}$ is less than threshold.

### 3.3 Extracting tree vertex

The segmented results using watershed image segmentation method is only the candidate tree objects, which are not necessarily the real tree objects. So we must identify if or not every classified tree objects is tree in the next step. And the NDVI produced from the spectral bands of Quickbird image is used to do this judgment process. The mean NDVI is used to classify the tree objects. The objects with mean NDVI value less than the threshold will be classified as non-tree, which will be deleted; otherwise, which should be classified as tree.

The highest point in an object is classified as tree vertex after this classification, these are used to estimate tree height. We use a 5×5 moving window to calculate the tree height if there more than on maximum value in one image object, and the mean DN value in this 5×5 window acts as the tree height in this case. If the maximum DN value is still equivalent in 5×5 window, we will extend window size until we can get the unique maximum DN value.

### 3.4 Estimating the height of suppressed trees using Weibull function

The height of suppressed trees can't be estimated using those existing methods [14], including the method described in this Section 3.3. So we should develop some method to do this. Weibull function is proposed by Weibull, who is the Swedish physicist, which is used widely in forestry measuring and research for the flexibility of its form. Weibull function is

LiDAR usually only can find the trees in the dominant, and the distribution of the underlying wood is difficult to obtain [14], to make use of the Weibull function for modeling and prediction and extract the height distribution of the lower wood. The Weibull function is proposed by the Swedish physicist Weibull, because is widely used in forestry research, it is the main method for the prediction of the height distribution and the diameter distribution [15]. The three-parameter Weibull function of the density function f(x) and distribution function F(x):
$$f(x) = \begin{cases} \frac{c}{b} \left( \frac{x-a}{b} \right)^{c-l} \times e \left( \frac{x-a}{b} \right)^{f} & x \ge a \\ 0 & x < a \end{cases} \qquad (9)$$

$$F(x) = 1 - e^{-\frac{x-a}{b}} \qquad (10)$$

a is a location parameter, b is a scaling function, c is a shape parameter.

In the application of using the Weibull function to simulate the trees high, a is the height of minimum spanning tree height, so that a becomes a set value .At this point, the three-parameter Weibull function is converted to two-parameter Weibull function used here are a set worthy of the two-parameter Weibull function. Because the determination of the laser radar, lidar accuracy combined with mutual shadowing between trees, and stacked cover, trees and undergrowth of wood is difficult to be identified, measured by the tree is not up to the minimum spanning tree, so only is the measured number of large trees, trees number is unknown, the corresponding can not be used to complete the Weibull function to the left truncated Weibull function to simulate the distribution of trees, get the value of b and c, then a complete weibull function simulate the distribution of the trees. Left truncated Weibull function of the density function f (x) as follows:

$$f(x) = \begin{cases} \frac{c}{b} \left(\frac{x}{b}\right)^{c-l} \times e\left(\left(\frac{t}{b}\right)^{c} - \left(\frac{x}{b}\right)^{c}\right) & x > t \\ 0 & x < t \end{cases} \qquad (11)$$

t is the measured height of the laser radar.

Weibull function parameter estimation There are several ways, the most common and accurate maximum likelihood estimates the parameters of the weibull function, 11 b and c, values, and then substituted into equation (10) to get a high degree of less than t probability F (t).

## 4 Results and analysis

Using the algorithms discussed in section 2 to segment LIDAR data of the two pilot areas, and the results as shown below.

<figure><img src="image_7.png" /><figcaption>Fig.6. Segmentation results ((a) segmentation result of test area one (b) segmentation result of test area two)</figcaption></figure>
We can be seen from Figure 6 that trees are basically of two pilot areas have been extracted, and the border in line with the actual situation, some spread between the artificial construction of the trees are extracted. However, due to the spectrum appear synonym spectrum with the spectrum of foreign bodies, and some buildings were classified as trees, some trees are not extracted.

NDVI can remove a large number of non-vegetation information. Can be seen from Figure 5, trees and buildings are well separated. But can be found in Figure 5 (a): the top right of the red marked man-made structures has also been judged as trees, the average NDVI of these areas is 0.4, but the bottom left of the green marked regional NDVI trees region NDVI only 0.36. Staggered addition of artificial buildings and trees, mixed distribution, but when the building height to a certain extent, the shadow of the building will be obscured to some trees, then trees NDVI values also fall, leading to an average of some trees image object NDVI dropped to less than the threshold, thus false positives for non-trees objects. Artificially built-up areas in the second image, the trees did not get to some of the buildings shadow, but in the forest part of the overall NDVI values are relatively high, so the shadow of the forest part of the less affected. In order to retain more information of the crown, set low NDVI, the first image of NDVI threshold is relatively high with 0.35, while the second image of 0.3.

Image segmentation and image objects of the merger after the completion of the judgment, the process to complete the Full Lambda the Schedule method, this method taking into account the spectral information of the Quickbird image object, area information and public side length of the adjacent area. Spectrum and size of the difference is proportional to the object with two images, side length and area with the adjacent edges of the two image objects is inversely proportional to the differences. This method is more reasonable to merge some of the region has a long common edge, but because of the high spectral similarity of the two image trees region, when the trees were over-segmentation, some of this should be a whole tree is divided into multiple regions ,but adjacent region area does not necessarily have similar spectral information, if it belongs to the edge of the area from the whole partition, it would be more similar to the edge of the area and other adjacent tree spectrum rather than it should be subordinate the tree itself. And because of the generally over-segmentation regions of various shapes, so the shape parameter and the height distribution of the information is difficult to join .So we use spectral information, but the merger of the two adjacent regions, the spectral information is too close to the merger, only two adjacent regional spectral information over a certain range when the merger. In this paper, spectral information is also used NDVI values.

Extracted from the tree height and the measured points were compared, the results are as follows:
<figure><img src="image_8.png" /><figcaption>Fig.7. tree height fitting results</figcaption></figure>

We can see in Figure 7 that using LiDAR measured tree height and measured
height has a high correlation, which is close to the 1:1 ratio. Measured tree height =
1.12 × the LiDAR tree height +0.087, R2 = 0.93.

Simulated trees using Weibull function, and predict the two images to detect the elevation distribution of the trees greater than 5m tree. 4057 and 1284, respectively, on the simulation results of the trees greater than 5m R20.96 and 0.90, predicted less than 5m of trees 26 and 12. As the first image is mainly forest areas, so using the tree height distribution measured in the first image as a test, set to detect minimum tree height is 15m, left-truncated Weibull function simulation to obtain between 5-15m the tree height distribution, the results show that the 3702 tree greater than 15m 5-15m of the tree 443, the distribution of more than 5m tree obtained by the Weibull function and to detect the distribution of greater than 5m tree compared to the R2 is 0.93. The above three Weibull function simulation through the Kolmogorov - Smirnov test can account for that using left truncation Weibull function can be effective simulation and forecasting detected and undetected trees.

# 5 Conclusions

This paper combined the use of LiDAR data with Quickbird images and constructed the model of finding tree height which take into account the height of the trees of lower wood. This method relative to the maximum search approach, be able to very good use of the spatial variation of the height distribution, rational use of trees and between trees existing gap in the watershed segmentation process, bifurcation identification of trees overlap or tree good effect; After the first watershed segmentation using watershed segmentation on a high level and some smooth crown can be split effectively, and using the Full Lambda Schedule algorithm can effectively merge the over-segmentation area, so not only guarantee the adequacy of the segmentation to reduce over-segmentation; Combination of LiDAR data and Quickbird images can effectively distinguish between forest and other surface features; Using the Weibull function to simulate the distribution of tree height and predict the distribution of the trees. This method of height information extraction can obtain high extraction efficiency in the region of dense forest and artificial buildings and trees interlaced distribution. Proposed in this study to consider the underlying
wood to a high degree of tree height to find the height model extraction process does not require too much prior knowledge, application scalability. Of course, the proposed algorithm there is need for further improvements:

1. In the process of tree height extraction, the threshold of NDVI and canopy size need to manually set and the threshold size would affect the results, the NDVI impact on the distinction between forest and other surface features, and crown size affects the crown dense region segmentation results.

2. In the process of canopy object merging, each region is relatively small and are independent, the use of the consistency of every region of the spectrum to the merger will bring some errors. If you can use the multi-echo or echo image of the laser radar, canopy and lower operating directly on point clouds, merged correctly.

## Acknowledgements

This research is funded by Agricultural Research Elitist Fund of Ministry of Agriculture of the People's Republic of China, Open-end Fund of Key Laboratory of Agri-informatics, Ministry of Agriculture, P. R. China [2012003], and National Key Technology R & D Program of the Ministry of Science and Technology [2012BAD20B0103]. We would like to thank the reviewers for valuable comments and constructive advice.

## References

1. Xiao, X. W. China's forest resource inventory. China Forestry Publishing House, Beijing (2005)

2. Andersen, H. E., Reutenbush, S. E., Mcgaughey, R. J. A rigorous assessment of tree height measurements obtained using airborne LiDAR and conventional field methods. Canadian Journal of Remote Sensing, 32(5):355--366 (2006)

3. David, R. S. Nancy, F. G. LiDAR measurement of sagebrush steppe vegetation heights. Remote Sensing of Environment. 102:135--145 (2006)

4. Husch, B, Miller, C, Beers, T. Forest mensuration 2th ed. Ronald Press Company, New York (1972)

5. Hans, O. Ø., Erik, N., Ole, M. B. Classifying species of individual trees by intensity and structure features derived from airborne laser scanner data[J]. Remote Sensing of Environment. 113: 1163--1174 (2009)

6. Liu, Q. W. Study on the Estimation Method of Forest Parameters Using Airbone LiDAR, Beijing (2009)

7. Persson, A., Holmgren, J., Soderman, U. Detecting and measuring individual trees using an airborne laser scanner. Photogrammetric Engineering & Remote Sensing, 68(9):925--932(2002)

8. Sorin, C. P., Randolph, H. W. Seeing the trees in the forest: using lidar and multispectral data fusion with local filtering and variable window size for estimating tree height. Photogrammetric Engineering & Remote Sensing, 70(5): 589--604 (2004)

9. Daniell, A. Z, David, L. E., George, C. Characterizing vertical forest structure using small-footprint airborne LiDAR. Remote Sensing of Environment, 87(2):171--182 (2003)

10. Tomas, B., Timothy, A. W., Rick, E. L. Detection and analysis of individual leaf-off tree crowns in small footprint, high sampling density lidar data from the eastern deciduous forest in North America. Remote Sensing of Environment, 85: 290--303(2003)
11. Solberg, S., Maesset, E., Bollandsas, O. M. Single tree segmentation using airborne laser scanner data in a structurally heterogeneous spruce forest. Photogrammetric Engineering & Remote Sensing, 72(12):1369--1378(2006)

12. Vincent, L, Soille, P. Watersheds in digital spaces: An efficient algorithm based on immersion simulations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(6):583--598(1991)

13. Robinson, D. J., Redding, N. J. Implementation of a fast algorithm for segmenting SAR imagery. Scientific and Technical Report, 1(1): 4--5(2002)

14. Maltamo, M., Eerikainen, K., Pitkanene, J. Estimation of timber volume and stem density based on scanning laser altimetry and expected tree size distribution functions. Remote Sensing of Environment, 90: 319--330(2004)

15. Meng, X. Y. A study of relation between D and H-Distribution by using the Weibull function. Journal of Beijing Forestry University, 1(10): 40--48 (1988)