Horvitz-Thompson-like estimators based on stochastic geometry for forest remote sensing

(1)

Dissertations in Forestry and Natural Sciences

KASPER KANSANEN

HORVITZ–THOMPSON-LIKE ESTIMATORS BASED ON STOCHASTIC GEOMETRY FOR FOREST REMOTE SENSING

PUBLICATIONS OF

THE UNIVERSITY OF EASTERN FINLAND

(2)

(3)

PUBLICATIONS OF THE UNIVERSITY OF EASTERN FINLAND DISSERTATIONS IN FORESTRY AND NATURAL SCIENCES

N:o 390

Kasper Kansanen

HORVITZ–THOMPSON-LIKE

ESTIMATORS BASED ON STOCHASTIC GEOMETRY FOR FOREST REMOTE

SENSING

ACADEMIC DISSERTATION

To be presented by the permission of the Faculty of Science and Forestry for public examination at the University of Eastern Finland, Joensuu, on November 6th, 2020,

at 12 o’clock.

University of Eastern Finland School of Computing

Joensuu 2020

(4)

Grano Oy Jyväskylä, 2020

Editors: Pertti Pasanen, Matti Vornanen, Jukka Tuomela, and Matti Tedre

Distribution:

University of Eastern Finland Library / Sales of publications julkaisumyynti@uef.fi

http://www.uef.fi/kirjasto

ISBN: 978-952-61-3492-5 (print) ISSNL: 1798-5668

ISSN: 1798-5668 ISBN: 978-952-61-3493-2 (pdf)

ISSNL: 1798-5668 ISSN: 1798-5676

(5)

Author’s address: University of Eastern Finland School of Computing

P.O.Box 111, 80101 JOENSUU, FINLAND email: kasper.kansanen@uef.fi

Supervisors: Professor Lauri Mehtätalo University of Eastern Finland School of Computing

P.O.Box 111, 80101 JOENSUU, FINLAND email: lauri.mehtatalo@uef.fi

Professor Matti Maltamo University of Eastern Finland School of Forest Sciences

P.O.Box 111, 80101 JOENSUU, FINLAND email: matti.maltamo@uef.fi

Associate Professor Aku Seppänen University of Eastern Finland Department of Applied Physics

P.O.Box 1627, 70211 KUOPIO, FINLAND email: aku.seppanen@uef.fi

Professor Jari Vauhkonen University of Helsinki

Department of Forest Sciences

P.O.Box 27, 00014 HELSINKI, FINLAND email: jari.vauhkonen@helsinki.fi Reviewers: Research Professor Erkki Tomppo

Aalto University

Department of Electronics and Nanoengineering P.O.Box 15500, 00076 AALTO, FINLAND email: erkki.tomppo@aalto.fi

Research Professor Johannes Breidenbach Norwegian Institute of Bioeconomy Research National Forest Inventory

P.O.Box 115, NO-1431 ÅS, NORWAY email: job@nibio.no

Opponent: Research Associate Thomas Opitz

French National Institute of Agronomic Research (INRAE) Biostatistics and Spatial Processes Unit

84914 AVIGNON, FRANCE email: thomas.opitz@inrae.fr

(6)

(7)

Kasper Kansanen

Horvitz–Thompson-like estimators based on stochastic geometry for forest remote sensing

Joensuu: University of Eastern Finland, 2020 Publications of the University of Eastern Finland Dissertations in Forestry and Natural Sciences

ABSTRACT

In this dissertation, Horvitz–Thompson-like estimation of forest characteristics from remote sensing data is studied. In forest remote sensing, the forest is measured with technical devices such as airborne or terrestrial laser scanners. These measurements produce three-dimensional point clouds, from which trees can be extracted with individual tree detection algorithms. These algorithms produce a sample of tree objects and tree attributes from the forest area that has been remotely sensed.

Possible attributes are, for example, height, crown diameter, and stem diameter at breast height. Usually all trees cannot be detected. There can be several reasons for this, but a major one is that some of the trees produce in some way nonvisible areas where other trees can be located and remain undetected. This partial detection produces underestimation of population totals and bias to mean and distribution estimates calculated directly from the attributes of the detected trees. In this thesis, individual tree detection is thought of as a sampling procedure where trees have different probabilities to be included in the sample based on their attributes.

If these probabilities can be approximated, Horvitz–Thompson-like estimators for forest characteristics of interest can be formed. Here we assume that the detection probabilities, or detectabilities, are related to the geometry of the pattern formed by the detected trees. We also assume that the forest is generated by a marked point process. These assumptions give us the tools to approach the calculation of detectabilities, and make the Horvitz–Thompson-like estimation possible. We study the estimation problem in two different remote sensing situations, namely airborne laser scanning and terrestrial laser scanning. We also study the problem of estimating attributes that have not been observed.

Universal Decimal Classification (UDC):519.246, 528.8, 630*58 Mathematical Subject Classification (MSC):62D05, 60G55

Library of Congress Subject Headings: Remote sensing; Stochastic geometry; Point processes; Mathematical models; Forests and forestry - Measurement

Yleinen suomalainen ontologia: kaukokartoitus, matemaattiset mallit, stokastiset proses- sit, metsänarviointi

Keywords: Point patterns, Forest inventory

(8)

ACKNOWLEDGEMENTS

The research for this study was carried out during 2014–2020 in the School of Com- puting at the University of Eastern Finland. This study was funded by the UEF spearhead projectMulti-scale Geospatial Analysis of Forest Ecosystems, the projectFour- dimensional Airborne Laser Scanningat the Department of Applied Physics, UEF, the project Towards Semi-supervised Characterization and Large-area Planning of Forest Re- sources Using Airborne Laser Scanning Data Acquired for Digital Elevation Modellingat the Department of Forest Sciences, University of Helsinki, the projectordSpat: Mod- els of Heterogeneity, Contextuality and Self-interaction in Ordered Spatial Point Patterns with Applications to Animal Movement and Forest Research, a consortium project between UEF and Natural Resources Institute Finland, the Finnish Cultural Founda- tion, North Karelia Regional fund, and The Doctoral Programme in Science, Tech- nology and Computing at UEF. I want to thank my supervisors Lauri Mehtätalo, Matti Maltamo, Aku Seppänen and Jari Vauhkonen and co-authors Petteri Packalen and Timo Lähivaara for their support and collaboration. I also want to thank Lauri, Matti, Jari and Aku for helping me to acquire funding from the aforementioned projects. I would like to thank the staff of the School of Computing for providing a pleasant work community. I would also like to give special thanks to Antti Pent- tinen, Juha Heikkinen, Mari Myllymäki and Mikko Kuronen for useful discussions and commentary on my work, and reviewers Erkki Tomppo and Johannes Breiden- bach for their commentary on this dissertation.

I thank my wife for her tremendous support during this process. This thesis is dedicated to her.

Helsinki, September 9th, 2020 Kasper Kansanen

(9)

LIST OF PUBLICATIONS

This thesis consists of the present review of the author’s work in estimation of forest characteristics from remote sensing data and the following selection of the author’s publications:

I K. Kansanen, J. Vauhkonen, T. Lähivaara, and L. Mehtätalo, "Stand density estimators based on individual tree detection and stochastic geometry,"Canadian Journal of Forest Research46(11), 1359–1366 (2016).

DOI: 10.1139/cjfr-2016-0181

II K. Kansanen, J. Vauhkonen, T. Lähivaara, A. Seppänen, M. Maltamo, and L.

Mehtätalo, "Estimating forest stand density and structure using Bayesian individual tree detection, stochastic geometry, and distribution matching,"ISPRS Journal of Photogrammetry and Remote Sensing152, 66–78 (2019).

DOI: 10.1016/j.isprsjprs.2019.04.007

III K. Kansanen, P. Packalen, M. Maltamo, and L. Mehtätalo, "Horvitz–Thompson- like estimation with distance-based detection probabilities for circular plot sampling of forests,"Biometrics(2020). Published online.

DOI: 10.1111/biom.13312

Throughout the overview, these papers will be referred to byI,IIandIII.

AUTHOR’S CONTRIBUTION

The publications selected for this dissertation are original, peer-reviewed research articles. The ideas for papersI-IIIoriginated from discussions between the author and co-authors, most significantly Lauri Mehtätalo in all cases, Jari Vauhkonen in case ofII, and Matti Maltamo in case ofIII. The author has carried out all numerical computations and necessary implementations of new methods in papers I-III. The author has carried out the analysis and discussion of results in papers I-IIIin collaboration with the co-authors. The author has written papersI-IIIas the main author, with significant contributions from the co-authors in all cases. The author has written the present review, which functions as an introduction to the methodology considered inI-III, with helpful commentary from the supervisors and reviewers.

(10)

(11)

1 INTRODUCTION

Forest inventory refers to the collection of forest information and data [1, p. xi].

Commonly collected information includes, for example, measurements of tree diameters and heights, and aggregated forest characteristics such as mean diameter, number of trees, and total volume of trees in some forest area. Examples of data used in forest inventory include field measurements, aerial photographs, satellite imagery, and laser scanner measurements, to name a few. Forest inventory can be seen as ecological and economic bookkeeping; having an up to date picture of the natural resources, and their ecological and economic value. Forest inventory can be used to support decisions relating to forest management, such as timing of thinnings or cuts. Many countries have national forest inventories [2, pp. 1–18].

Information from a national forest inventory can be used on one hand to formulate regulations and recommendations for forest management to ensure sufficient sup- ply of resources for forest industries, and on the other hand to plan forest protection actions on national level and to negotiate and monitor global agreements such as the the Kyoto Protocol [3, pp. 1–16].

Forest remote sensing refers to the collection of forest information without direct measurements. The examples of forest inventory data mentioned above, excluding the field measurements, belong to this category. Remote sensing supported forest inventory also includes the endeavour of deriving information on tree and/or forest characteristics from this remote sensing data. What makes forest remote sensing interesting is the ability to collect data from larger areas in shorter amount of time when compared to traditional inventory methods. Airborne and space-borne methods make what is called ”wall-to-wall” measurements possible; this means that an entire forest area can be covered with measurements, whereas field measurements in contrast are usually gathered only from smaller field sample plots located in different parts of the forest area. Airborne and space-borne methods also make it possible to collect data from inaccessible areas where movement on the ground is hard, not allowed (e.g. military zones), or too dangerous (e.g. areas that have been covered with anti-personnel mines).

Two main methodologies for derivation of forest characteristics from remote sensing data are the area-based approach (ABA) and individual tree detection (ITD).

For a review of these methodologies in relation to airborne laser scanning, see [4].

In ABA, summary statistics of remote sensing data that correlate with aggregate forest characteristics are calculated over areas where field measurements and remote sensing data are available. Then, either a parametric model – e.g. a linear regression model – is fitted between the statistics and values of characteristics and used to predict values of the characteristic in areas where field measurements are not available, or a nonparametric method such ask-nearest neighbours is used for the prediction. On the other hand, ITD tries to find tree objects from the remote sensing data, producing tree-level predictions. Naturally, the tree-level predictions can be aggregated to form predictions of forest characteristics, for example mean height. It should be noted that ITD cannot be used with all remote sensing data;

low resolution data, such as very sparse laser scanning measurements or satellite

(14)

images, do not contain enough information to detect individual trees. It should also be noted that the terms ABA and ITD are usually used in conjunction with airborne laser scanning. However, there is no conceptual difference between using summary statistics of laser return height measurements to explain aggregate forest characteristics, and using summary statistics derived from e.g. aerial photographs to explain those forest characteristics (see e.g. [5]). Similarly, there is no conceptual difference between searching individual trees from laser point clouds and searching those trees e.g. from aerial photographs (see e.g. [6]). Hence, the terms ABA and ITD are used in this work, regardless of the type of the remote sensing data, to describe approaches that conceptually fit to these descriptions.

This thesis focuses on ITD from airborne and terrestrial laser scanning. ITD can be used to derive tree-level information, and some forest characteristics such as number of trees per unit area – the stem density – seem to be more directly connected to the collection of detected tree objects than the summary statistics used in ABA, at least at a conceptual level. As an illustrative example, let us consider stem density estimation more carefully. In ABA, one could, for example, fit a (Poisson) regression between the number of trees and some covariates derived from remote sensing data. This leads to a model selection problem, where such covariates that best explain the variation in stem densities between different areas have to be chosen. There is, of course, no guarantee that any of the possible covariates are related to stem density, although the number of possible covariates from which to choose is usually very high. On the other hand, ITD will produce a stem density estimate directly – the number of detected trees.

Although ITD is a very promising methodology, it does have its problems. The problem that this thesis focuses on is that not all trees can usually be detected. The reason behind these nondetections is that the detection of a tree is affected by the other trees. It is easy to imagine, for example, that when the forest is observed from above, the tallest trees with the largest tree crowns can be spotted very easily, but smaller trees located under the large crowns are not visible. This can lead to underestimation in some forest characteristics such as number of trees, and overes- timation in others such as mean height, because the small trees are missing from the calculation. Another error in ITD that can occur is the algorithm interpreting a group of close, small trees as one big tree. This could be seen as one of the trees being detected correctly, some attributes estimated erroneously, and that tree caus- ing the nondetection of the others. An error to the other direction is also possible:

the algorithm interpreting one large tree as several small ones. Commonly, the errors where an existing tree has not been detected are calledomissionerrors and the errors where nonexisting trees are detected ascommissionerrors. For example, the comparison of different ITD methods in [7] showed omission errors ranging roughly from 25 to 60 per cent of the existing trees, and commission errors ranging roughly from 25 to 55 per cent of the detected trees, leading to total number of detected trees ranging roughly from 50 to 140 per cent of the existing trees. The systematic errors produced by ITD can lead to erroneous management decisions, which can have significant economic consequences [8]. This thesis concentrates on the correction of errors in forest characteristic calculations produced by the omission of smaller trees located under larger tree crowns, and does not present methods for correcting the other type of omission error or commission errors.

ITD can be seen as a sampling algorithm. Some trees are detected, i.e. included in a sample, and others are not. It is natural to assume that some properties of a tree, e.g. its size and location with respect to other trees, have an effect on being

(15)

detected or not. More precisely, there exists a probability for a tree with certain attributes to be detected based on all the locations it could have been located at, and this probability is affected not only by the attributes of the tree in question, but also by the other trees. Horvitz and Thompson [9] provided an estimator for sampling situations where the probabilities of objects to be included in a sample vary. They also showed that the Horvitz–Thompson estimator is unbiased. There exist also formulas for the variance of the estimator and approximate confidence intervals [10, p. 70]. Horvitz–Thompson estimator assumes that the sampling probabilities are known beforehand, which is of course not the case with ITD. However, Horvitz–

Thompson-like estimators, where the sampling probabilities are approximated from the observed data, are commonly used for example in distance sampling, where total sizes of animal populations are estimated from the number of detected animals [11].

Horvitz–Thompson-like estimators, and the problem of approximating the detection probabilities, ordetectabilities, from ITD data are in the heart of this thesis.

The randomness related to the detectabilities, as mentioned above, arises from the thought that the trees are located randomly in the forest. They form a pattern that is the result of complicated natural or human-induced process, which can be modelled with a stochastic process that governs the locations and attributes of trees. These processes are known as marked point processes, and the study of them and all random patterns with certain geometrical properties is known as stochastic geometry[12]. The tools provided by stochastic geometry are the basis for deriving detectabilities in this thesis.

ITD does not necessarily give estimates of all tree attributes of interest. For example, ITD from airborne remote sensing data most likely gives information on tree heights and crown sizes, but it does not necessarily contain information on the diameters of tree stems at breast height. In these cases, the Horvitz–Thompson-like estimation has to be enhanced with transformations from the detected tree attributes to the attributes of interest. In this thesis, we discuss the use of distribution matching with mixed-effects models to extend the Horvitz–Thompson-like estimation to this problem.

The purpose of this thesis is to present Horvitz–Thompson-like estimators with detection probabilities based on stochastic geometry for ITD-based forest inventory.

The main objective of this work is to present methods for mitigating the problems of ITD mentioned above, hence improving the operational usability of ITD as a tool of forest inventory. This thesis functions as a concise introduction to the concepts and methodologies used in I-III. Although the introduction is brief, it is broader and more unified than what has been written in any single publication included in this thesis. The intention behind this thesis is that a reader without any knowledge in forest inventory, remote sensing or stochastic geometry can get a good overview of the subject matter as a whole. First, basic concepts relating to forest inventory, forest remote sensing and stochastic geometry are reviewed. After this, methodology relating to the Horvitz–Thompson-like estimators and approximation of detection probabilities based on stochastic geometry are presented, after which methods for estimating unobserved forest characteristics are introduced. The overview is concluded with a summary of results ofI-IIIand discussion on possible future research topics related to the current work.

(16)

(17)

2 INTRODUCTION TO FOREST INVENTORY

Information relating to practical forest inventory can be found in [1], especially inventory based on field measurements. This section is based on [1, pp. 3–134].

Suppose there exists a forest area with N trees, indexed by i = 1, . . . ,N, with attributes Y_i^(k), e.g. stem diameter, stem volume or biomass. The purpose of forest inventory is to find population totals

Y^(k)=

∑

N i=1

Y_i^(k), (2.1)

means

Y¯^(k)= ^Y^(k)

N , (2.2)

ratios

R_k,l= ^Y

(k)

Y^(l), (2.3)

and other statistics such as variances, covariances and correlations. Examples of population totals, means and ratios are number of trees, mean diameter of tree stems at breast height and the proportion of certain tree species, respectively. Other examples of forest characteristics that are generally of interest are basal area

BA=

∑

N i=1

π

DBH_i 2

2

, (2.4)

where DBH is the stem diameter at breast height – this is a population total of the cross-sectional area of trees at breast height – and quadratic mean diameter

QMD=

vu ut

∑

^N

i=1

DBH²_i

N , (2.5)

the square root of the squared mean stem diameter. In addition, the whole distribution of DBH is a valuable characteristic [13], preferably by tree species [14]. It should be noted that we use Nhere to denote the total number of trees. However, typically in literature relating to forestry Nis scaled to trees per hectare.

In practice, the forest area can be quite large and have a very large number of trees N. Hence, it is not possible or economically viable to measure every tree, and moreover,Nis not usually known. From these factors it follows that the forest inventory is based on sampling and estimation. The sampling strategy is based on distributing several sampling points on the forest area and measuring trees around these sampling points. This measurement produces afield plot. A commonly used type of a field plot is thefixed-radius circular sample plot, which can be established by running a string of length equal to the desired plot radius from the sampling point, the plot centre point, outwards [15]. Sampling strategy based on circular sample

(18)

plots is known as circular plot sampling. In addition to circular plots, square and rectangular field plots are also commonly used. The distribution of sample points can followsimple random sampling, where the sampling points are located completely randomly to the forest area, orclustered random sampling, where clusters of sampling points are located randomly to the forest area, or systematic sampling, where the sampling points follow a grid with a fixed distance between sampling points and the location of the starting point of the grid and the directionality of the grid are randomized. The sampling strategy can also combine the aforementioned strategies:

for example, clusters of sampling points can be located in the inventory area based on a systematic grid.

Forest inventories can be categorized toone-phaseandtwo-phasesampling schemes.

In a one-phase scheme only the field sampling is performed. In a two-phase scheme the first phase consists of collecting auxiliary information from a larger number of sampling points, and the second phase consists of field sampling from a subsam- ple of the points sampled in the first phase. Remote sensing supported inventories could be included into this two-phase scheme, where the first phase of auxiliary information collection refers to the remote sensing data collection and the second phase corresponds to the gathering of field data for possible calibration of a model that connects the remote sensing data to the underlying forest characteristics.

There are two conceptually different approaches to estimating forest characteristics from samples, commonly called design-basedinference and model-based inference [16]. In design-based inference, it is assumed that the population is fixed, and the sample is random. In other words, randomness related to the inclusion of certain trees in the sample comes from the random locations of the sampling points, which depends on the design of the sampling scheme. The inference does not depend on the distributions of the forest characteristics in this case. Contrastingly, model-based inference assumes that the locations of trees, their number Nand attributesY_i^(k)are random, and hence the randomness related to the sampling of certain trees follows from the random process governing the forest. This process is represented by a model, and the validity of the inference depends on how well the chosen model represents reality.

Sometimes different terminology is used to categorize these two frameworks for inference. Hansen et al. [17] definesprobability-sampling designto encompass a sampling plan that produces known probabilities > 0 of inclusion in the sample for every member of the population, and estimators for which the validity does not depend on an assumed model, at least when the sample size is large enough. They definemodel-dependent sampling designas any design where either the sampling plan or estimators are chosen because they have desirable properties under an assumed population model, e.g. unbiased estimators under model assumptions. They consider model-based designs to be different from model-dependent ones; for example a design where the sampling procedure or estimators are based on a model, but the design coincides otherwise with probability sampling in the sense that the validity of the results does not depend on the model, is not model-dependent. It is useful to understand that conceptually the design-based framework of [16] and probability-sampling design of [17] correspond to each other, as do the model-based framework of [16] and the model-dependent design of [17], although the terminology differs. Additionally, Gregoire [16] remarks that the semantic distinction of model-based and model-dependent approaches, made by Hansen et al. [17], is not widely used and can lead to confusion. In this thesis we use the terms model-based

(19)

and model-dependent as synonyms to describe approaches where the validity of model assumptions affects the validity of estimation results and inferences.

The model-based framework can have advantages when compared to the design- based framework, for example useful inferences from small samples and possibility to search for a best estimator, leading to smaller sampling errors than what would be possible with other estimators [17]. These advantages of course depend on the underlying model corresponding to reality. In this thesis, a model-dependent approach is adopted, but not necessarily because of the previously mentioned reasons.

Rather, the sampling problem that is considered, the sampling of trees by ITD from a tree population, is different from the typical sampling situation described above.

What is considered here is the estimation of forest characteristics on a connected forest area – this could range from a plot to the full inventory area – when only some of the trees on that area are sampled, i.e. detected, by ITD. In this case there is no possibility to design a sampling plan that would produce known detection probabilities for all of the trees – the detection is related to the random process of forest and the trees affecting the detection of each other, not a predefined sampling plan. Hence, a model-dependent approach must be used. This leads to the use of Horvitz–Thompson-like estimators, whereas in the design-based framework the Horvitz–Thompson estimator [9] could be used due to the known detection probabilities.

(20)

(21)

3 FOREST REMOTE SENSING

This chapter gives a short introduction to two data collection methods of forest remote sensing most relevant for the work reviewed here, airborne laser scanning (ALS) and terrestrial laser scanning (TLS). Both methods are based onLight Detection and Ranging(LiDAR) technology, where a scanner emits light, the light reflects back from an object as anechoand is read by the scanner system. The scanner measures the time that it takes for the emitted light to return back to the scanner and based on the known speed of light produces a measurement of the distance from the scanner to the object. Lastly, the chapter gives a short introduction to Individual Tree Detection (ITD), a methodology for interpreting forest remote sensing data.

3.1 AIRBORNE LASER SCANNING

This section is heavily based on [4], which is a good introduction to the forestry applications of ALS and a good source of references for further reading. In ALS, a LiDAR scanner is mounted on an aircraft. The position and orientation of the scanner are recorded using Global Positioning System (GPS) and Inertial Navigation System (INS). The light emitted from the scanner mostly hits the forest canopy and the ground. Based on the known position and orientation of the scanner and the distance of the object from the scanner, the location and height of every hit can be calculated. Hence, ALS produces 3D surface measurements of the forest canopy and the ground.

There are two types of LiDAR scanners: discrete return and full waveform. The difference between these two systems is how they handle the recording of the data produced by an emitted laser pulse. The forest canopy can be seen as a porous material. Hence, it is unlikely that all of the energy contained in a laser pulse would be returned from the same height. It is more likely that different amounts of the laser pulse penetrate the tree crowns to different depths, even to the ground level below the canopy. The discrete return systems only record one or a few return echoes for every emitted laser pulse, the first one being the highest measurement of the canopy, the last one being close to the ground level. Full waveform systems, on the other hand, digitize the full sequence of returned energy and produce a functional representation of the return energy, as a function of height. The discrete return systems are more commonly used in operational remote sensing supported forest inventories, and most useful are the first return heights, representing the surface of the canopy, and the last return heights that can be used to form a digital terrain map that can further be used to remove the effect of varying terrain elevation from the laser data. Figure 3.1 presents an example of this kind of normalized ALS point cloud.

ALS data collection has some parameters that can affect the quality of the data.

Parameters related to the actual aircraft are the flying altitude and speed. Parame- ters related to the scanning device are beam divergence, pulse repetition frequency, scanning pattern, and opening angle. Beam divergence describes how the diameter of the laser beam increases as the distance from the scanner increases. Together

(22)

Figure 3.1: Airborne laser scanning data, and individual tree detection data derived from it.

with the flying altitude it determines the footprintof the laser, the diameter of the laser beam on the ground, basically the area from which a laser pulse produces measurements. Small footprints – 1 metre or less – are often preferred in forestry applications [4, p. 5]. Pulse repetition frequency refers to the amount of pulses emitted by the scanner per time unit. Opening angle refers to the angle at which the scanner scans from side to side, producing a corridor of scanned area around the flight path, the width of which depends on the opening angle and flying altitude.

High flying altitude and large opening angle makes it possible to scan large areas with a single flight path. However, a large opening angle will produce measurements more from one side of the trees – somewhat like a profile – near the edges of the aforementioned corridor. All of the parameters mentioned earlier will affect the density of the ALS point cloud, the number of points per unit area. Wang et al. [7] concluded that point density of 2 points·^m⁻²can provide satisfactory results if ITD is used to only measure dominant trees. However, they also anticipated that larger point densities would improve the characterization of 3D canopy structure.

For example, the ALS data used inIandIIhas a nominal sampling density of 11.9 pulses·^m⁻²^.

In forestry, the strength of ALS is that the 3D structure of the forest canopy can be accurately measured and characterized. In addition, very large spatial areas can be covered quite easily, wall-to-wall.

3.2 TERRESTRIAL LASER SCANNING

This thesis considers TLS as a method of circular plot sampling, where the LiDAR scanner remains stationary during the scan. It is also possible to mount the scanner on a vehicle and scan along some path, see e.g. [18]. A good review of TLS in the context of forestry is [19], which is used as a source for this section. For circular plot sampling, the scanner is placed on a sampling point, standing on a tripod.

The scanner rotates 360 degrees in the horizontal plane, while it scans vertically

(23)

depending on an opening angle. This produces a 3D point cloud that accurately characterizes the forest structure below the canopy around the scanning point – at least the parts that are not occluded by trees or other obstacles close to the scanner.

If only one scan from a circular plot is performed, this setup is called a single- scan setup. However, it is also possible to try and avoid the occlusion problems by scanning the same plot from different points outside the plot in addition to the centre point scan. This is called amulti-scansetup. Naturally, the multi-scan setup increases field measurement time, and more work is needed to combine the data from the different scans.

TLS has similar scanner parameters as ALS, and the scanner can record discretely one return, the first and the last return, some additional returns between the first and the last, or full waveform. In addition to LiDAR based on the return time of the laser pulse, instruments based on the phase-shift of the laser are used. These systems record only the first return, have very high point-densities and fast data- acquisition speeds, and are suitable for measurement of objects less than 100 metres away from the scanner. The instruments based on the return time are slower, but can be used at larger distances, and offer more recorded returns per laser pulse. In general the pulse densities produced by TLS instruments are high to accommodate the detection of tree objects from the point clouds; for example, Litkey et al. [20]

utilized a scanner that produced measurements 6 mm apart at a distance of 10 metres from the scanner.

3.3 INDIVIDUAL TREE DETECTION

ITD, also sometimes called individual tree crown (ITC) detection or delineation in the airborne remote sensing context, is a methodology that tries to find individual tree objects from remote sensing data and estimate tree-level attributes for the detected trees. Several methods have been developed both for ALS and TLS data.

In ALS-based ITD the tree objects are usually tree crowns, whereas in TLS-based ITD the tree objects are usually tree stems. It should be noted that usually in the literature ITD refers to the ALS-based methods, and the term is not used in conjunction with TLS. It is the choice of the author to unite the task of finding tree objects from different sources of remote sensing data under one umbrella, ITD, as the term clearly characterizes the methodology and its independence of the data source. Let us first concentrate on the ALS case.

One quite simple image analysis based method, presented by [21] and [22], takes a canopy height model formed by interpolating the laser data – more specifically, the first and only return heights – and finds local height maxima through low- pass filtering. These local maxima are interpreted as locations of trees. After this, the canopy height model is segmented into tree crown objects around the local maxima via watershed segmentation. An estimate of crown diameter is achieved by taking the maximum of the diameters of the corresponding segment in four cardinal directions passing the crown location. This method has been used inIIfor benchmarking.

Quite different approach is presented by Lähivaara et al. [23], the method used inIand IIfor ITD. This method models a tree crown as a parametric, rotationally symmetric geometrical object. The objects are parametrized by the horizontal coor- dinates of the tree crown centre point, crown radius, crown height, the lower limit of the living crown, and the crown shape. It is assumed that the ALS data is produced

(24)

by a generative model that consists of a collection of these geometrical objects and error terms related to the porosity of the canopy and the related penetration of the laser pulses, and the fact that tree crowns are not rotationally symmetric, smooth geometrical objects. Based on this model, the problem can be studied in the framework of Bayesian inverse problems [24]. Field measurements and allometric models of [25] can be used to form a prior distribution for the model parameters. The fi- nal product is the maximum a posteriori (MAP) estimate of the model parameters, which represents a collection of geometrical objects representing the tree crowns that can be detected from the ALS data. An example of ITD with this method is shown in Figure 3.1.

Another 3D method is presented in [26], where the laser returns from a tree crown are modelled directly as a probability distribution and the returns from an area as a mixture of the tree-specific distributions within the area. Other methods, mostly working in similar vein to [21], search for local height maxima but additionally use allometric knowledge [27], a priori knowledge of expected crown size or tree density [28,29], or multilayered approach [30] to improve the detection of trees.

Recently, ITD methods based on deep learning have also emerged, see e.g. [31]. It should be noted that ITD algorithms are quite active area of research, and there are many algorithms in addition to the examples given here. For example, Zhen et al. [32] reviewed 212 articles of research literature related to ITD from 1990 to 2015, 92 of which concentrated on algorithm development.

For TLS, most of the ITD methods utilize geometrical models. For example, Lovell et al. [33] first takes a slice of the TLS data at certain height and classifies measurements as belonging to tree stems if their apparent reflectance exceeds a threshold value. Angular point intensities are used to infer stem locations. Circu- lar stem models relating reflectance values to the stem diameter are then applied.

Liang et al. [34] first classifies TLS data points as either belonging to tree stems or not based on the properties of local covariance matrices of point locations. Then, the stem points are divided into clusters, each cluster containing points from one stem, based on distances between the points. Cylindrical models are fitted to the point clusters. Raumonen et al. [35] uses what they callquantitative structure models, which represent trees as hierarchical collections of cylinders or other shapes. In addition to stem modelling, these models are also able to represent the branches.

Ravaglia et al. [36] uses the Hough transform and growing open active contours, also called snakes, to model trees as tubular shapes consisting of series of 3D circles.

InIII, which considers estimation in the TLS case, no actual TLS data was used, and hence none of the aforementioned ITD methods were used. This study was purely a simulation study, exploring the behaviour of Horvitz–Thompson-like estimation in different spatial patterns and comparing the estimator to two other previously proposed estimators.

(25)

4 STOCHASTIC GEOMETRY

This chapter gives a short introduction to the tools of stochastic geometry that are central to the calculation of detection probabilities presented in I–III. First, point processes, models of random locations and randomly located objects, are considered.

Second,mathematical morphology, a methodology of transforming geometrical structures to obtain more information on those structures, is presented.

4.1 POINT PROCESSES

This section is based on [12] and [37]. Point processes are models for irregular point patterns. Very general definitions could be formulated, but we constrain ourselves and only consider point processes defined in a compact subset of the planeS⊂^R²^. Apoint processis a random variableXthat generates values fromN(S), the set of all finite discrete subsets ofS. One realisation ofXtakes the formx ={x₁,x2, . . . ,x_N}^, where x_i are points in S, x_i 6= x_j for all i 6= j, and the number of points N can be random. Some processes can be characterized via a probability density function describing the probabilities of different point configurations. However, for some point processes the density function cannot be given in analytical form, and in some cases it is more convenient to describe the process via some other properties. A point process is called stationary if the properties of the process do not depend on the location of the window W ⊂ S where the process is observed. Expressed differently, the distribution of the process is translation-invariant. An example of a nonstationary process would be a process where the number of points increases as we move from one side of S to the other. A point process is called isotropicif it is rotation-invariant; in other words, the window W can be rotated an arbitrary amount, and the observed properties of the process in the rotated and unrotated window should be the same. In this thesis it is always assumed that the underlying forest process is stationary and isotropic.

An important point process model is the (homogeneous)Poisson process. Some- times patterns generated by the Poisson process are said to exhibit complete spatial randomness(CSR). The Poisson process is characterized by two properties: first, the number of points located in a bounded Borel setB⊂Sfollows the Poisson distribution with expected valueλ|B|^{, where}|^.|is the area operator andλ>0 is called the intensity parameter of the process, and second, the numbers of points inkdisjoint Borel sets are independent of each other. These properties lead to the points of the process being independently located inS, so that there are no interactions between them. The Poisson process is especially important as a null model when testing if a pattern exhibits interactions between the points.

A special case of a point process is the single random point. For a random point x that is uniformly distributed in S, the probability that it hits a subset A is P(x ∈ A) = ^|A|_|_S_|. In other words, the probability is directly proportional to the area of the subset. The single random point has a strong connection to the Poisson process. Poisson process can be thought of as a collection of uniformly distributed random points.

(26)

In the literatureλis used generally to denote the intensity or global point density of a stationary point process, the expected number of points generated by the process per unit area. Mathematically, this can be expressed as

λ= ^E[N(S)]

|S| ^, ^(4.1)

where N(S)is the number of points in setS.

4.1.1 Characterizing point patterns and processes

Point patterns and processes are usually roughly categorized to three groups based on the interactions they exhibit. First category is the CSR, with no interactions.

Second category is regularpatterns. Regular patterns have less variation in numbers of points located in disjoint subsets of S than CSR and there are less of small inter-point distances than in a CSR pattern. An example of a regular process is a hard core process, where all the points are at least distance rapart from each other.

Third category isclusteredpatterns. These have more variation in numbers of points located in disjoint subsets ofSthan CSR and there are more points with small inter- point distances than in a CSR pattern. An example of a clustered pattern is a point pattern generated by first generating parent points via a Poisson process and then generating the actual points around these points in some local neighbourhoods, e.g.

limiting the distance from a child point to a parent to be less or equal to some upper limitr.

Point patterns can exhibit different spatial characteristics at different scales. For example, tree centre points will always be some minimum distance apart, deter- mined by stem size, but the trees could still be growing in groups. The point pattern formed by the tree centre points would hence exhibit regularity at short distances and clustering at longer distances. Due to this scale-dependence point patterns are often characterized by functional summary statistics, statistics that are functions of distance. Two commonly used statistics are (Ripley’s) K-function and (Besag’s)L- function.

TheK-function compares the local point density in the neighbourhood of a typical point of the process to the global point density:

K(r) = ^E[N(B(o,r)\ {o})]

λ , r≥^0, ^(4.2)

where

B(x,r) ={y∈^R²^:||x−y|| ≤r}^. ^(4.3) is the x-centred disc of radius r. Here it is assumed that the typical point of the point process is located on the origino, which is possible due to stationarity. For the Poisson process

K(r) =πr², (4.4)

because from the definition of the Poisson process it follows that

E[N(B(o,r)\ {o})] =λπr². (4.5) Clustered processes have more points aggregated near a typical point than the Pois- son process, which leads to larger local point densities thanλand

K(r)>πr². (4.6)

(27)

On the other hand, regular processes have more isolated points than the Poisson process and

K(r)<πr². (4.7)

TheK-function can be estimated from an observed point pattern by calculating the mean number of points in the neighbourhoods of all the points, divided by an estimate of λ. However, the limited size of observation window W can produce bias to the estimated function – the full neighbourhoods of points near the window boundary are not observed. Hence, edge-effect corrections are necessary. The size of Walso affects the maximum distancerat which theK-function should be estimated – for example, in a circular window W the maximum r commonly used is half of the radius of the window. Estimators of theK-function with several edge-effect corrections and sensible defaults of maximum distances can be found in the package spatstat[38], implemented in R [39].

TheL-function is defined via theK-function as L(r) =

rK(r)

π , r≥^0. ^(4.8)

Hence, in the Poisson process case L(r) = r, and for clustered processes L(r) > r and regular processes L(r) < r. This normalization makes visual inspection of an estimated function with respect to a theoretical function easier. The L-function is estimated through the estimator ofK-function. Examples of point patterns and their K-functions andL-functions are presented in Figure 4.1.

A measure of deviation from CSR can be formulated via theL-function. First, the estimated bL(r) is calculated for distances r up to some maximum distance R.

Then, the distance at which the absolute difference between the estimate and the theoretical function is the largest is located:

r^∗=arg max

r∈[0,R]|r−^bL(r)|^. ^(4.9)

After this, the signed difference at this distance is used as the deviation measure:

r^∗−^bL(r^∗). (4.10)

This measure is close to zero if the pattern exhibits CSR. Positive values indicate stronger signs of regularity than clustering, and negative values indicate stronger signs of clustering than regularity. This measure was used in III to measure de- viations of point patterns – and collections of patterns via means of estimated L- functions – from CSR.

It should be noted that there exist nonfunctional summary statistics for characterizing spatial point patterns, such as the Clark-Evans index [40]. This index is the ratio between the mean nearest-neighbour distance in the observed pattern and the expected nearest-neighbour distance of a Poisson process with the same intensity as the observed pattern. The Clark-Evans index was used inIand IIto measure the deviation of patterns from CSR.

4.1.2 Marked point processes

An important extension of point processes aremarked point processeswith realisations {(xi,mi)}_i=1^N ^{. Here,}xiare locations, just like in point processes, but now there is an

(28)

Cluster CSR Regular

point patternKfunctionLfunction

function theoretical estimated

Figure 4.1: Examples of point patterns exhibiting three different characteristics:

clustering, complete spatial randomness and regularity, and theK- andL-functions calculated from the point patterns. The number of points and the size of the window is the same in all cases.

(29)

additional markmi – possibly generated by a random variable – connected to each of the point locations. Marked point processes also have concepts of stationarity and isotropy, the distributional invariance covering in this case also the marks. The mark can be any attribute of interest, for example diameter, height, volume, or biomass of a tree. There can be several different marks connected to one point; for example, a point could be marked simultaneously with all of the aforementioned examples of marks. The mark can also be a geometrical object: a commonly considered geometrical mark in the literature is the disc B(x,r). Marked point processes with geometrical marks are commonly called germ-grain modelsin the literature. Germ- grain models are usually used as models of random sets by taking a union over the geometrical marks. For example, ifx_iare the locations of points as above and every point has as a mark a disc with radiusr_i, then

[N i=1

B(xi,ri) (4.11)

is the random set generated by that process. One of the most studied germ-grain models is theBoolean model, where the underlying point process is the Poisson process and marks are independent identically distributed random compact sets independent of the Poisson process, for example discs with random radii. The Boolean model has been used e.g. to model the spatial pattern of heather [41]; several other applications have been listed in [12].

Germ-grain models can be characterized by several statistics emerging from their connection to sets. In the planar case, standard measures of sets are for example the surface area covered by the set, the boundary length of the set, and the number of connected components of the set. The germ-grain model produces a random set;

hence, these measures are also random. Due to this, the expected values of the set measures are a natural way to characterize the germ-grain models. However, very often in practice only one realisation of the model is observed in an observation win- dowW⊂S. In some special cases, such as the Boolean model, formulas connecting expected values of the set measures to model parameters are known. Generally, however, these connections are not known. The two aforementioned points can make the fitting of germ-grain models a challenging task. In this thesis the fitting of germ-grain models is not of interest. The germ-grain models are simply a tool upon which the presented estimators are built. Nevertheless, it is good to keep in mind that deriving information relating to an underlying process from a single realisation has its difficulties. This relates to the concept ofergodicity: if a process is ergodic, a single, large enough sample can be analysed to get statistically meaningful results.

A sufficient condition for ergodicity is mixing, which in a sense means that distant parts of a (marked) point process are independent.

4.2 MATHEMATICAL MORPHOLOGY

Mathematical morphology is a methodology for transforming geometrical structures to other geometrical structures with certain central operators. Very general theory can be formed; see e.g. [42, pp. 13–35]. For our purposes, we can concentrate on morphology of subsets ofR²; see e.g. [12, pp. 6–8]. Usually the morphological transformations are combined with measurements of the transformed set. For example, a set is transformed several times, and every time, the area, boundary length, or the number of connected components of the transformed set is calculated. This

(30)

dilation erosion

Figure 4.2: Examples of the morphological transformations erosion and dilation performed on a planar set with the black disc.

procedure could be compared to extending a linear regression model with transformed covariates to gain additional information on the process being modelled.

As an example, Arns et al. [43] used measures of transformed sets to characterize different types of disordered materials.

The morphological operators are built upon structuring elements K ⊂ ^R²^{. Dif-} ferent shapes of K lead to different types of transformations and hence different interpretations of measurements made from the transformed sets. To simplify the situation we only consider the origin centred disc B(o,r) = {x : ||x|| ≤ r} ^{as a} structuring element. Theerosionof setA⊂^R²with the discB(o,r)is the set

A⊖B(o,r) ={x∈ A:B(x,r)⊂A}^. ^(4.12) In other words, the erosion of a set A with B(o,r) is the set of those points in A where the centre point of a disc of radius rcan be located so that the disc is fully contained inA. Thedilationof setAwith the discB(o,r)is the set

A⊕B(o,r) ={x∈^R²^:B(x,r)∩A6=^∅}^. ^(4.13) In other words, the dilation ofAwithB(o,r)contains all the points where the centre point of a disc of radius r can be located in such a way that the disc comes into contact withA. Erosion and dilation of a planar set with a disc have been presented in Figure 4.2.

(31)

5 HORVITZ–THOMPSON-LIKE ESTIMATORS

Horvitz and Thompson [9] presented an estimator for population totals when sampling without replacement from a finite population. The estimator is a sum of the sampled attributes of interest, over which the estimation of a total is wanted, weighted with the probabilities of sampling the subjects that were sampled. For the Horvitz–Thompson estimator, the sampling probabilities are known beforehand.

However, when ITD is used for tree detection, and in a sense for sampling the underlying tree population, the probability of detecting a certain tree, the detectability, is not known. The probabilities must be estimated somehow. To make the difference clear, estimators that use estimated detectabilities are called Horvitz–Thompson-like e.g. in distance sampling of animal population sizes [11]. These estimators were also calledHorvitz–Thompson typeinI. There are different approaches to calculating the detection probabilities, but our approach is based on interpreting forest as a realization of a marked point process, the assumption that the detection of a tree is affected by other trees, and certain assumptions relating to mathematical morphology. Other approaches have been studied e.g. by [44] and [45].

LetX = {(xi,mi,si)}_i=1^N be a realization of a marked point process onW ⊂ ^R² with attributes of interestmi and secondary markssi that affect the detectability of pointi. We assume that the point locations and marks do not contain measurement errors. We are interested in estimating a population totalτ=∑_i=1^N mi. The Horvitz–

Thompson-like estimator of this total can be written as

b τ=

∑

N i=1

miDi

p(s_i)^, ^(5.1)

where Di is an indicator variable of detection – 1 if pointi has been detected and zero otherwise – and p is the detectability that depends on s_i. During estimation, the sum is calculated over the detected points, and hence the marks or detectabilities of undetected points are not needed.

To calculate the detectabilities, it is assumed that the detection of points has a sequentialnature. The sequence emerges from the simple idea that in the application cases considered in this thesis, a tree that is closer to the remote sensing device has an ability to affect the detection of a tree that is further away – for terrestrial remote sensing the closeness relates to the actual distance of the tree from the device, and for airborne remote sensing it relates to the height of the tree. More generally, this means that point i affects the detection of point j only if i < j. Nothing affects the detection of the first point, and hence it always has detectability 1. For point j>1 the pointsi =1, . . . ,j−1 produce areasAi⊂^R²with properties that depend on the location xi and mark si that affect the detectability of point j. To be more precise, a point with a mark sj would not be detected if xj ∈ G(^S_i=1^j⁻¹Ai), where Gis a morphological transformation that depends on sj. If we assume that xj is a single random point that is uniformly distributed inW the probability of not being

(32)

detected can be written as

P



x_j∈ G





j−1 [ i=1

A_i



^\W



=

G_S_j−1

i=1Ai

_T W

|W| ^, ^(5.2)

and hence the detectability is

p(sj) =1−P



xj∈G





j−1 [ i=1

Ai



^\W



. (5.3)

It should be noted that the detectabilitiesp(s_j)that we have formulated are con- ditional probabilities, conditioned on the points i < j, in one sense the observed

”history” of the process. Assuming that the observationsi< jare fixed and x_j is a single uniformly distributed random point, the detectabilityp(s_j)is not an estimate, but a known fixed probability.

It is also possible to condition the detectability, if it is known that xj must belong to a certain subset of W, by changingW to that subset in the aforementioned formulas. The subset can be of lower dimensionality thanW, for example a line or a circle – a situation encountered inIII. In this case|^.|measures length. It is useful to think Xas a point process in a larger window Sthat containsW, so that all the points xi that lie outsideW but have such an area Ai that comes into contact with W can be taken into account during detectability calculations, as these points can affect the detectability of points inW. Depending on the application and available data the detectability can depend on all points i = 1, . . . ,j−1, if it is possible to derive Ai also for undetected points from the data, or only those points that have been detected. In the latter case Ai = ∅ for the undetected trees in Equations 5.2 and 5.3.

Horvitz and Thompson [9] showed that the Horvitz–Thompson estimator is unbiased, i.e. E[bτ] =τ. Similarly, Horvitz–Thompson-like estimator is unbiased if the assumptions related to the detection process and underlying marked point process hold; an example proof is presented inIIIfor data generated by the Poisson process.

For the Horvitz–Thompson estimator, formula for the variance is known [10, p. 70]:

var(_bτ) =

∑

N i=1

1−pi

p_i

m²_i +

∑

N

i=1

∑

j6=i

pij−pipj

p_ip_j

!

mimj. (5.4)

This theoretical formula depends on all of the members of the population N and cannot actually be used for variance calculation. However, an unbiased estimator of the variance, calculated with the detected members of population, is also known [10, p. 70]:

varc(τ_b) =

∑

idetected

1 p²_i − _p¹

i

!

m²_i +2

∑

idetected

∑

jdetected,j>i

1 p_ip_j − _p¹

ij

!

mimj, (5.5)

where pi = p(si) and pij is the probability to include bothi and j in the sample.

This estimator of variance can be used with the Horvitz–Thompson-like estimator and should hold when the aforementioned model assumptions hold.

(33)

An important special case of the Horvitz–Thompson-like estimator is the estimator of the population size, i.e. the number of pointsN, which is given by settingb mi=1 for alli:

Nb =

∑

N i=1

Di

p(s_i) =

∑

idetected

1

p(s_i)^. ^(5.6)

For the detected pointsithe estimatorNb can be divided into componentsNbi= _p(s¹_i₎, so that Nb = ∑Nbi. Nbi can be interpreted as an estimator of the number of points that are similar to point i. This leads to an estimate of mark distribution, and an estimate of the cumulative distribution function can be given for s– assuming the mark is one-dimensional – as

Fb(s) = ¹ Nb

∑

s_i≤s

Nb_i, (5.7)

where the indexing is again over the detected points. This formulation can also be extended to give distribution estimates over other observed marks of interest by replacings_iwith the mark in question in Equation 5.7. TheNb_ialso make it possible to estimate the mean of a mark as a weighted mean

1 Nb

∑

idetected

Nbimi. (5.8)

5.1 DERIVING DETECTABILITIES FROM AIRBORNE REMOTE SENS- ING DATA

InI, the problem of estimating stem density – that is, the total number of trees in a forest area – from ALS based ITD was considered. For stem density estimation mi =1. Although the algorithm of Lähivaara et al. [23] which was used for ITD in the paper produces 3D tree objects, the situation was simplified by considering the problem only in two dimensions. This led to si being the maximum crown radius of treei, andAi= B(xi,si). It was assumed that only the trees with larger value of s_icould affect the detection of treejwiths_j <s_i. This led to the sequence of points being ordered according tos_i, so thats1>s2>. . .>sN. It was assumed that a tree would not be detected if a certain proportion of its crown was covered by the larger trees. This proportion was controlled by introducing a tuning parameter α ≥ ^{0 to} the morphological transformation, and the detectability became

p(sj) =1− h_S_j₋₁

i=1A_i

⊖B(o,αs_j)ⁱ^TW

|W| ^. ^(5.9)

The interpretation of parameter α is as follows. If α = 1, it is assumed that a tree will not be detected if the crown of the tree is fully covered by the larger trees. If α = 0 – meaning no erosion is performed – it is assumed that a tree will not be detected if the centre point of the crown is covered by the larger trees. A value of αbetween 0 and 1 means that at least α per cent of the crown radius has to be covered for detection to fail. The limiting case here is the centre point being covered. Naturally, this construction could be extended to allow also trees with crown centre points outside of the set of larger tree crowns to be left undetected by

(34)

Figure 5.1: An example of the morphological aspects of detectability calculation.

On the left, an example forest area with detected trees. Detectability is calculated for the tree marked as a black disc. On the center, the data are limited to contain only trees that are larger than the black one. On the right, an appropriate morphological transformation is performed, in this case erosion withα=1, producing a set where the tree could have been located in such a way that it would have been completely covered by the larger tree crowns.

replacing the erosion with a dilation. However, this was not considered inIorII. In this case the effect of undetected trees to the calculation of the detectability was not considered, andAi =∅for undetected trees. The derivation of the detectability has been illustrated in Figure 5.1.

5.2 DERIVING DETECTABILITIES FROM TERRESTRIAL REMOTE SENS- ING DATA

In III, the problem of estimating forest characteristics of interest from ITD based on single-scan TLS was considered. It was assumed that tree stems closer to the scanner could affect the detection of trees further away. To simplify the problem, it was assumed that the stems were discs with diameters equal to stem diameters at breast height d_i, and all of the effects on detection were related to these discs.

Heres_i = (r_i,d_i), wherer_i is the distance of the treeicentre point from the origin, where the scanner is placed. The assumption of trees closer to the scanner affecting the detection of trees further away leads to the points being ordered according to r_i−^d₂ⁱ, the closest distance from the scanner to the stem disc, so thatr1−^d₂¹ <. . .<

rN− ^d₂^N. Furthermore, this ordering quite naturally leads to considering the trees as random points uniformly distributed on a circle of radiusr_i,S(r_i), not the whole window W. In this case the A_i are formed by the stem disc and the area behind it, the subset of plane that is between the tangent lines of the stem disc that travel through the origin (see Figures 1 and 2 in III). The formula for the detectability became

Horvitz-Thompson-like estimators based on stochastic geometry for forest remote sensing

Dissertations in Forestry and Natural Sciences

KASPER KANSANEN

HORVITZ–THOMPSON-LIKE ESTIMATORS BASED ON STOCHASTIC GEOMETRY FOR FOREST REMOTE SENSING

Kasper Kansanen

HORVITZ–THOMPSON-LIKE

ESTIMATORS BASED ON STOCHASTIC GEOMETRY FOR FOREST REMOTE

SENSING

TABLE OF CONTENTS

1 INTRODUCTION

2 INTRODUCTION TO FOREST INVENTORY

∑

∑

∑

3 FOREST REMOTE SENSING

4 STOCHASTIC GEOMETRY

dilation erosion

5 HORVITZ–THOMPSON-LIKE ESTIMATORS

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑

∑