Characterizing Tree Spatial Distribution Patterns Using Discrete Aerial Lidar Data

: Tree spatial distribution patterns such as random, regular, and clustered play a crucial role in numerical simulations of carbon and water cycles and energy exchanges between forest ecosystems and the atmosphere. An efficient approach is needed to characterize tree spatial distribution patterns quantitatively. This study aims to employ increasingly available aerial laser scanning (ALS) data to capture individual tree locations and further characterize their spatial distribution patterns at the landscape or regional levels. First, we use the pair correlation function to identify the categories (i.e., random, regular, and clustered) of tree spatial distribution patterns, and then determine the unknown parameters of statistical models used for approximating each tree spatial distribution pattern using ALS-based metrics. After applying the proposed method in both natural and urban forest sites, our results show that ALS-based tree crown radii can capture 58 %-77 % (p < 0.001) variations of visual-based measurements depending on forest types and densities. The root mean squared errors (RMSEs) of ALS-based tree locations increase from 1.46 m to 2.51 m as the forest densities increasing. The Poisson, soft-core, and hybrid-Gibbs point processes are determined as the optimal models to approximate random, regular, and clustered tree spatial distribution patterns, respectively. This work provides a solid foundation for improving the simulation accuracy of forest canopy bidirectional reflectance distribution function (BRDF) and further obtain a better understanding of the processes of carbon and water cycles of forest ecosystems.


Introduction
Tree spatial distribution patterns, as an essential characteristic of forest structure, play a vital role in forest succession, regeneration, growth, and understory development [1][2][3][4][5][6][7]. It thoroughly explains the coexistence relationship and community structure of plant species based on past spatial distribution and ecosystem processes [8,9]. Accurate and quantitative representations of tree spatial distribution patterns are the critical foundation in studying forest radiative transfer processes [10], which further affect the simulation accuracy of carbon and water cycles using process-based ecological models. For example, tree spatial distribution patterns are essential inputs to geometric optical models to establish the relationships between forest structural parameters and forest canopy directional reflectance [11,12]. Thus, better representations of tree spatial distribution patterns will improve the accuracy of estimating the proportions of visible forest overstory and ground components in a forest stand, which further simulate the forest canopy bidirectional reflectance distribution function (BRDF) and retrieve forest structural parameters such as leaf area index (LAI). The spatial distribution pattern of trees is illustrated based on the relative tree locations in a horizontal plane [13], which included random, regular, clustered, or their combinations spatial distribution [14][15][16]. The same spatial coverage of a forest stand might exhibit different tree spatial distribution patterns under various sizes of observation windows [17].
During the past decades, researchers have developed some commonly used statistical models to approximate tree spatial distribution patterns based on tree location and crown size. For example, the Poisson point process is a well-accepted model used to approximate randomly distributed tree crowns at the landscape level [12]. The regular distributions are usually represented using the Gibbs point process [13,18], Markov point process [16], soft-core point process model [19], and Hypergeometric model [20]. In the case of clustered distribution patterns, the Neyman-A process [11], Marked area interaction point process [21], and Thomas process [22] are the most commonly used statistical models in practice. However, it is still challenging to determine the unknown parameters of various statistical models for a given specific forest stand in the landscape level due to the following reasons: (1) field data collection. The labor-intensive and time-consuming nature of fieldwork makes it challenging to obtain the tree locations and crown size information at the landscape level efficiently and accurately. Tree location can be represented by either stem location [23], treetop point [24], or crown center point [25,26]. Various ground-based instruments have also been used to collect tree stem locations at the plot level such as terrestrial laser scanning [27], total station theodolite [28], or rangefinder [29]. (2) tree crown competition: the competitive interactions [25] between tree crowns, characterized by the crown size information, are indispensable to determine tree spatial distribution patterns since tree crowns cannot approach each other infinitely close. However, few studies have paid attention to characterize competitive interactions between trees in acquiring limited resources (i.e., light, water, and nutrients) in the clustered tree spatial distribution patterns [11,17]. Thus, a generic mathematical model is needed to account for this effect based on crown information [30]. The hybrid-Gibbs model holds the ability to explain different distribution patterns at various spatial scales [31,32]. It can be used to approximate clustered tree spatial distribution patterns considering the competitive interaction between tree crowns.
The optical remotely sensed imageries are traditionally used to determine tree locations at the landscape level based on the local maxima filtering or photogrammetry-based visual interpretation methods [33][34][35]. However, due to the complexities of three-dimensional (3-D) canopy structure and topographical variations, it is challenging to map the crown size accurately at the landscape level using only optical remotely sensed data [36], especially in heterogeneously mixed forests. With the advancement of remote sensing technology, the detailed 3-D structural information implicitly contained within light detection and ranging (lidar) data [37,38] makes it possible to capture directly tree spatial distribution patterns [39]. Aerial laser scanning (ALS) data have been successfully used to estimate tree heights, treetop locations, and crown sizes using various algorithms [40][41][42][43]. For example, Packalén et al. [39] identified types of tree spatial distribution patterns using ALS-based treetop locations without traditional statistical models. Then, a powerful and flexible ALS-based approach should be developed to quantify various tree spatial distribution patterns. The goals of this study were: (1) Develop an ALS-based method to characterize tree spatial distribution patterns at the landscape level quantitatively; (2) Characterize the competitive effects among tree crowns using ALS-based metrics; (3) Investigate the factors impacting ALS-based tree crown segmentation algorithms and the impacts of observation window sizes on tree spatial distribution patterns.
In this study, we first derived the tree locations and crown size information for two distinct urban heterogeneous and natural homogeneous forest stands based on the ALS data. Then, the point pair correlation function was used to identify the tree spatial distribution patterns. For each category, we approximated them using appropriate statistical models and determined their unknown parameters based on the ALS-based metrics. Finally, ALS-determined statistical models were validated using the Monte-Carlo simulation method. The flowchart of the current study was shown as follows ( Figure 1): Figure 1. The flowchart of characterizing tree spatial distribution patterns using aerial laser scanning (ALS) data.

Study Sites
Two study sites are selected in this study to test the performances of the proposed ALS-based approach in characterizing tree spatial distribution patterns at the landscape level. The first one locates in the Washington Park Arboretum (WPA), which is a well-managed urban forest. It is at the south of the University of Washington campus in Seattle, Washington, U.S.A. The dominant tree species of the WPA are douglas fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla), western redcedar (Thuja plicata), big-leaf maple (Acer macrophyllum), monkey puzzle (Araucaria araucana), southern magnolia (Magnolia grandiflora), and new mexican locust (Robinia neomexicana). The slope throughout the WPA site is less than 15 % (Figure 2). The second study site is at the Panther Creek (PC) watershed, which is a natural but still managed forest on the east side of the coastal mountain of Oregon, USA. It is located 57 km southeast of Portland, Oregon, USA at 45° 18′ N, 123° 21′ W. The elevation of Panther Creek watershed ranges from 100 m to 700 m. The douglas fir (Pseudotsuga menziesii), western hemlock (Tsuga heterophylla), western redcedar (Thuja plicata), and big-leaf maple (Acer macrophyllum) are dominant tree species. The spatial patterns at these two locations, due to the management practices and natural processes, are expected to be distinguishable and capture a wide range of tree spatial distribution patterns.

Datasets
 Aerial laser scanning data (ALS) The airborne lidar data for the WPA site were acquired using a Riegl LMS-Q560 laser scanner (Riegl, Co. Ltd. Horn, Austria). The coverage was collected at an average elevation of 279 m and flight height was about 310 m on average. The scan angle of the laser scanner ranged from −30° to 30° from nadir. This ALS flight setting resulted in a discrete point dataset containing five returns per pulse and more than ten pulses per square meter.
In the PC site, the laser sensor Leica ALS80 (Leica Geosystems AG, Heerbrugg, Switzerland) mounted in an aircraft with a frequency of 394.8 kHz emitted a single laser pulse and view field of 30° at the flight height of 1400 m. Finally, it yielded a 3-D lidar data with an average density of greater than eight pulses per square meter over terrestrial surfaces. The average point densities in WPA and PC sites are 26 points / m² and 15 points / m², respectively.  Validation data To validate tree crown segmentation result and applicability of ALS-determined statistical models, we set up four squared forest plots in the WPA site, including a high-density coniferous plot (i.e., WPA-HC), a medium-density coniferous plot (i.e., WPA-MC), and low-density coniferous and mixed forest plots (i.e., WPA-LC and WPA-LM) with a side length of 100 m. Furthermore, two squared forest plots with the side length of 100 m in the PC site were also selected, including one high-density coniferous plot (i.e., PC-HC) and a medium-density coniferous forest plot (i.e., PC-MC). In the six forest plots, we visually identified treetop points and obtained total tree numbers and their locations from ALS point cloud data. By doing this, forest plots were divided into high-(> 100 trees), medium-(≥ 80 trees and ≤ 100 trees), and low-densities (< 80 trees) groups. The crown diameters were measured in crossed directions (North−South and East−West) for each tree crown in ALS-based forest data. Moreover, we identified forest types visually using aerial photographs with a spatial resolution of 0.3 m. More detailed information about the forest plots is summarized in Table 1.  Modeled data To illustrate the significant effects of tree spatial distribution patterns on estimating the probabilities of viewing the ground surface in a forest scene, we modeled six square forest plots with the side length of 100 m and identical individual tree ALS point cloud data (Table 1). We first simulated three different tree spatial distribution patterns (i.e., Random, Clustered, and Regular) in the R programming environment [44], then explored the probabilities of viewing the ground surface area from various directions. Moreover, each spatial distribution pattern was parameterized using two different sets of parameters.

ALS-Based Tree Crown Segmentation
We obtained the locations and crown sizes of each tree using the previously published ALS segmentation algorithm [41]. The spacing threshold and search radius are two primary parameters in the algorithm. To overcome the over-or under segmentation issue, we developed a "point cloud slicing-based tree verification" method to achieve a better segmentation result. In this study, we used the treetop point as tree location due to the inability to capture stem location using aerial remotely sensed data.
For the point cloud data of each segmented tree, we extracted six sliced planes with a thickness of 0.2 m in two different ways horizontally and vertically, respectively. The six slicing planes located at the 15 %, 30 %, 45 %, 60 %, 75 %, and 90 % heights (i.e., horizontal slicing) or distances along the X-or Y-axis (i.e., vertical slicing) of the segmented individual tree points. We computed the standard deviations of the X-and Y-coordinates (i.e., horizontal slicing), or X-and Z-, or Y-and Z-coordinates (i.e., vertical slicing) of the six center points from six sliced plane using the Equation 1: We determined the threshold as 0.4 m based on the distance deviations between treetop and bottom points in the case of a horizontal slicing way. As for the vertical slicing way, the threshold was set as 0.8 m, the one half of the minimum tree crown size by assuming tree crowns are symmetrical. By finishing the iteration process, we obtained the final results of the treetop location, crown diameter, and tree height at the individual tree level.

Tree Spatial Distribution Pattern Identification
The point pair correlation function (i.e., g(r) function) is a second-order summary statistic function determining the types of tree spatial distribution patterns (i.e., random, clustered, or regular) for a given point set [17,45]. The g(r) function is calculated based on the number of point pairs for a given arbitrary point within a distance of r [45]. It is proportional to the first-order derivative of Ripley's K(r) function [46] concerning r: It contains the same statistical information as the K(r) function but the non-cumulative effect of K(r) function [47]. (1) Random distribution: if g(r) is constant and equal to one (i.e., g(r) = 1), it means that the location of any point is completely independent of the other point locations. It belongs to a completely random process (CSR) [48]; (2) Clustered distribution: if values of g(r) are greater than one (i.e., g(r) > 1), this suggests a clustering process caused by the aggregation effects of points; (3) Regular distribution: if values of g(r) are smaller than one (i.e., g(r) < 1). It indicates that the spatial distribution patterns of points tend to a regular pattern. It should be noted that there is a "virtual aggregation" phenomenon [17] in a forested area with patches of non-forest (i.e., human activities or deforestation) land at the landscape or region scales. In this case, the g(r) function suggests a clustered tree spatial distribution pattern across all spatial scales [17], and it is necessary to remove the nonforest areas to determine tree spatial distribution patterns.
To analyze the visual-and ALS-based g(r) function curves effectively, we computed the normalized g(r)norm function as where g(r)norm is the normalized value of g(r) function ranging from −1 to 1; Gmax(r) and Gmin(r) represent the simulation envelopes obtained as the 5 th highest and lowest values after 999 simulations using the Monte Carlo simulation method [49]. This normalization procedure allows us comparing different envelope lines more efficiently.  [50] among tree crowns quantitatively, we proposed a new individual tree crown model. In this model, a tree crown was represented using two annulus circles with treetop points as the center point (i.e., Pi). The radii of two circles were HCi and Ri, respectively ( Figure 3a). The "hard-core" (i.e., inner circle) was a circular area of individual tree crown where no parts of other trees existed, while the "buffer region" (i.e., outer annulus circle) area might have overlapped parts between tree crowns. In the current study, we used the hard-core radius to describe the competitive effects among tree crowns: the larger hard-core radius, the more substantial competitive effect.

ALS-Based Statistical Model Parameters Determination
According to relative spatial locations among tree crowns, the interactions among tree crowns are classified into three categories: (a) separated tree crowns-two tree crowns (i.e., Tree-1 and Tree-2) with the treetop points P1 and P2, hard-core radii and buffer circle as HC1 and HC2, R1 and R2, respectively ( Figure 3b). These are apart from each other when the distance between two treetop points is larger than the sum of radii of two tree crowns (i.e.,‖P1-P2‖> R1+ R2); (b) tangent tree crowns-two tree crowns may be adjacent to each other (i.e., they share a tangent line) when the distance between two treetop points is equal to the sum of the radii of two buffer circles (i.e.,‖P1-P2‖= R1+ R2, (Figure 3c); (c) overlapped tree crowns-two or more overlapped tree crowns within the tree crown "buffer regions" based on the geometric model of tree crowns (Figure 3d). For example, Tree-1 is overlapped in the buffer region by Tree-3 and Tree-4 whose treetop points and buffer regions radii are P3, P4, and R3, R4, respectively. The distances between every two treetop points are smaller than the sum of any two crown radii (i.e., ‖P1 -P4‖ < R1 + R4 and‖P1 -P3‖ < R1 + R3). , hard-core distance (i.e., HC1), buffer region radius (i.e., R1) (inset a), and three possible relative location relationships between tree crowns including two separated tree crowns (inset b), two tangent tree crowns (inset c), and three overlapped tree crowns (inset d) in practice. R13 and R14 represent the overlap distance between trees.

Hard-core radius determination
Determining the "hard-core" radius (HCi) is a prerequisite to quantify the competitive effects among tree crowns using the hybrid-Gibbs model. In the current study, we obtained the hard-core radius of a forest plot using the proposed individual tree crown model. For an individual tree, it may be overlapped with multiple tree crowns. In the case of two tree crowns, we computed the hard-core radius of the i-th tree crown overlapped with the j-th tree crown by assuming that the hard-core radius of a tree crown was proportional to their crown sizes as where HCij is the estimated hard-core radius when the i-th tree crown overlapped with the j-th tree crown; Ri and Rj are the radii of the i-th and j-th tree crowns; Rij is the line segment between the treetop points of the i-th and j-th tree crowns falling in the "buffer" regions and computed as ( ) If the target tree (i.e., the i-th tree) was overlapped by two other tree crowns (i.e., j-th and k-th tree crowns), we computed the hard-core radius of the i-th tree (i.e., HCi) by averaging the estimated hardcore radii (i.e., HCij and HCik). For example, we first computed the overlapping distances (i.e., R13 and R14) based on Equation 5, and obtained the estimated hard-core radii (i.e., HC13 and HC14) for Tree-1 concerning Tree-3 and Tree-4 using Equation 4 ( Figure 3d). The hard-core radius of Tree-1 was computed as the average values of HC13 and HC14. We used the averaged hard-core radii for all trees as the hard-core for the whole forest plot.

Random distribution:
We approximated a random tree distribution pattern using the simplest and well-accepted Poisson model [11,12,51] and computed as where is the average tree number per unit area; P(y) is the probability of finding y trees in the sample area. The unknown parameter could be determined using the ALS-based metrics.

Regular distribution:
We used the "soft-core point process model" to approximate a regular tree spatial distribution pattern. Compared to hard-core point processes, the interaction of points in the soft-core point process is smoothly decreasing with distance increasing. It is an essential point process for biology because there is rarely an abruptly ending for the interactions between points [52]. The "soft-core" model of a regular distribution point set x = {x1, x2, … xn} could be computed as where β is a contributing factor of each point to the global point process; σ controls the distance at which the interactions occur and is a positive real number. The larger values, the stronger interactions; κ indicates the strength of interactions between a point pair, and it is dimensionless variables ranging from 0 to 1; ‖xi − xj‖ is the distance between the two points. β is computed using the maximum likelihood estimator method [53], κ and σ are estimated using the maximum profile pseudolikelihood algorithm [54] using the ALS-based tree locations.

Clustered distribution:
The Neyman type-A and hybrid Gibbs models are the two most commonly used statistical models to approximate clustered tree spatial distribution patterns [11,31]. For comparison purposes, we have analyzed these two models in the current study. In a Neyman type-A model, there are "parent" (i.e., center points of every group) and "offspring" points (i.e., all non-parent points). If there are i trees in the given j groups within the quadrat, the probability of having i trees in the given j groups can be determined using these double-Poisson point processes and computed as where m1 is the mean number of groups in the quadrat, and m2 is the mean size of the groups. In the process of model realizing, m1 is approximated by the density of group centers ( N), and m2 is represented by group radius (rN) together with the mean number of points per group (mN). The minimum contrast estimation method [55] could be used to determine the unknown parameters ( N, rN, mN) of the Neyman-A model using ALS-based tree locations. However, the Neyman-A model only treats tree crowns as an individual point without considering the crown size and their interactions. Hybrid Gibbs model can approximate spatial distribution patterns under different scales using various pairwise interaction functions. Then we built a hybrid Gibbs model by combining the "hardcore" and "Geyer" point processes [56] to better characterize clustered tree spatial distribution patterns by incorporating the competitive effects among tree crowns (i.e., hard-core radius). The probability density of a point sets based on the newly constructed hybrid Gibbs model can be computed as where r is the distance between a point pair; H(r) is the probability density of the "hard-core" point process model; x is an unordered point set x ={x1, x2, … xn}; n is the total number of points in the point set (n ≥ 0); s is the saturation threshold to limit the number of distinct unordered points that are closer than r units apart; t(xi, x\xi, rg) is the number of neighbors of xi in point set x; rg represent the interaction ranges, and γ controls the strength of the repulsion (i.e., γ < 1) and attraction (i.e., γ > 1) interactions between points in the process. The parameter β controls the density of a point process; for a given γ, the density of point process increases as β increasing. The value of β equals the density of the Poisson process when γ = 1. The values of s and γ can be estimated using the maximum profile pseudolikelihood algorithm based on ALS-based tree locations [54], β can be obtained using the maximum likelihood estimator method [53]. We determined the hard-core radius using the method described in section 2.5.1 based on ALS-based tree locations and crown sizes. The values of H(r) and P(x) equal 0 when r is smaller than 2×HC, suggesting that no point is generated. H(r) equals 1 when r is greater than the 2×HC. In this case, the hybrid Gibbs model will be equivalent to the Geyer model. We can rewrite Equation 9 as By doing this, we can determine the unknown parameters of statistical models used for approximating various tree spatial distribution patterns quantitatively.

ALS-Determined Tree Spatial Distribution Models Validation
To test the suitability of the ALS-determined statistical model for a given forest plot, we compared the 95 % confidence intervals simulated from the ALS-determined statistical model with the "reference g(r) function curve." The reference g(r) function curve was produced based on the tree locations determined through visual interpretation. We generated the envelope lines through repeated random sampling (i.e., typically 999 times) based on the ALS-determined statistical models using the Monte-Carlo simulation method [49]. If the reference g(r) function curve fell into the envelope lines generated using a specific ALS-determined statistical model completely, the corresponding statistical model will be identified as a qualified model to approximate tree spatial distribution patterns.

Sensitivity Analysis
We compared the ALS-based tree locations and crown sizes with the ones obtained from the visual interpretation method by changing the spacing thresholds and search radii from 1 m to 5 m with the step of 1 m for the natural forest plots. Moreover, to investigate the effects of overlapping distances on the accuracy of tree crown segmentation, we conducted a numeric simulation experiment by changing the overlapping distances from 1 m to 4 m with the step of 1m in three different tree species configurations (i.e., conifer + conifer, broadleaf + broadleaf, conifer + broadleaf). In the forest plot WPA-HC, by changing the radii of hard-core from 1.4 m to 2.0 m with the step of 0.1 m, we explored the effects of hard-core on tree spatial distribution patterns.
Moreover, to illustrate the importance of tree spatial distribution patterns on simulating forest canopy directional reflectance using the existing method [57], we computed the probabilities of viewing the ground surface from various directions in principle plane. The azimuthal and zenith angles of sun were set as 90°and 30°, while the viewing azimuthal (VA) and viewing zenith angle (VZ) of observations directions were set as 90° and 0°, 90° and 30°, 90° and 60°, 270° and 30°, 270° and 60°, respectively.

ALS-Based Tree Locations And Crown Sizes
After applying the ALS-based tree segmentation algorithm, raw point cloud data were grouped into individual trees represented by different colors at in both WPA and PC sites.   By comparing the visual-and ALS-based tree crown radii ( Figure 6), we found that the values of the correlation coefficient increased from 0.61 to 0.77 as forest densities decreased in urban forest plots WPA-HC, WPA-MC, and WPA-LC (Figure 6a-c). Moreover, the mixed forest types were lower the value of R 2 in the low-density forest plots WPA-LC (i.e., 0.77) and WPA-LM (i.e., 0.66) (Figure 6c,  d). As for the natural forest stand, an increasing trend in R 2 was observed for the natural forest plots PC-HC (i.e., 0.58) and PC-MC (i.e., 0.71). Although, there were similar R 2 values between the ALSand visual-based tree crown radii in forest plots WPA-MC (i.e., 0.72) and PC-MC (i.e., 0.71), apparent RMSE differences were observed in WPA-MC (i.e., 1.39 m) and PC-MC (i.e., 0.85 m). Different forest vertical structures might explain this discrepancy in these two forest sites. Since the WPA site is a mixed urban forest with well-managed practices, it has an apparent adjacent configuration of overand under-stories with different tree species and crown shapes. However, in the PC natural coniferous forest, there are relatively fewer understory crowns. As shown in Table 2, it was found that in the urban forest canopies, the RMSE of ALS-based tree locations increased from 1.46 m to 2.51 m as the forest densities increasing. For the low-density forest plot, both the coniferous (i.e., WPA-LC) and mixed (i.e., WPA-LM) forest plots have small RMSEs of 1.46 m and 1.47 m, respectively. In the natural forest canopies, both high-and medium-density forest plots achieved small RMSE values for ALS-based tree locations due to the well-separated tree crowns in a natural coniferous forest environment. The mixed tree species and complicated vertical structure in the WPA site might increase the overlapping between tree crowns and further introduce errors into the estimation of tree locations.

ALS-Based Tree Spatial Distribution Patterns
Based on the ALS-and visual-based tree locations, we obtained two normalized g(r) functions for six forest plots in both WPA and PC sites (Figure 7). Overall, the ALS-based g(r)norm function curves have similar variation patterns with the visual-based ones for all forest plots. In forest plot WPA-HC (Figure 7a), the curves of both ALS-and visual-based g(r)norm functions were out of the normalized envelope lines. The ALS-based values of g(r)norm function exhibited different patterns (i.e., 1 < g(r)norm when 2.94 m < r < 6.65 m or g(r)norm < 1 when r < 1.81 m) at different ranges, indicating a mixed point spatial distribution at different spatial scales. However, the g(r)norm function of visual-based tree locations was smaller than 1 when r < 1.28 m, and larger than 1 when 1.96 m < r < 3.82 m. The tree spatial distribution pattern in WPA-HC was classified as a clustered distribution incorporated with competitive effects between tree crowns. In this case, a single scale statistical model, such as a Neyman-A model, is not enough for comprehensively characterizing the tree spatial distribution patterns. Thus, a generic and more sophisticated model, such as the hybrid Gibbs model, should be implemented instead.
Regarding the forest plot WPA-MC (Figure 7b), only part of g(r)norm function stayed within the normalized 95 % confidence intervals when the distance between pair points r was beyond 6.5 m. The visual-and ALS-based g(r)norm functions were out of the range of normalized envelope lines when distance r was less than 5.7 m and 6.5 m, respectively, indicating that the tree spatial distribution pattern was regular. Therefore, the soft-core point process model would be the optimum one to approximate the tree spatial distribution patterns in this case. As for forest plots WPA-LC and WPA-LM (Figure 7c, d), the ALS-and visual-based g(r)norm functions fell well within the normalized 95 % envelope lines (i.e., g(r)norm = 1 and −1), suggesting that the trees within these forest plots exhibited a random distribution pattern. Accordingly, the Poisson model should be used to approximate its tree spatial distribution patterns.
In the PC site, the forest plots PC-HC ( Figure 7e) and PC-MC (Figure 7f) have similar g(r)norm function curves with the forest plots WPA-MC and WPA-LM, respectively. In the forest plot PC-HC, the visual-and ALS-based g(r)norm functions were out of the normalized envelope lines when the r < 4.19 m and 5.02 m, respectively. The forest plots PC-HC and PC-MC were classified into regular and random distribution patterns. Thus, the soft-core point process model and the Poisson model were chosen to approximate the tree spatial distribution patterns. For the whole study area of both WPA and PC sites, the mix of forest and non-forest made the ALS-based g(r) function curves deviated from the 95 % confidence interval lines of the null model (i.e., CSR) (Figure 8a, c). Moreover, the extent of ALS-based g(r) function curves deviating from envelope lines in the WPA site was more significant than in the PC site because the WPA site is a well-managed urban forest with more non-forest area coverages. After masking out the non-forest land cover, we obtained a reasonable ALS-based g(r) function curves that partially fell into the envelope lines (Figure 8b, d). It was shown that the ALS-based g(r) were smaller than 1 when the r < 6.2 m and r < 7.1 m for WPA and PC sites, respectively. Therefore, the tree spatial distribution patterns in the WPA and PC sites were considered as regular distribution based on the criteria described in section 2.4, and the soft-core point process model was used to model the tree spatial distribution pattern. Figure 8. Comparisons between the g(r) function curves of the whole (i.e., forest + non-forest) (insets a and c) and forest only (insets b and d) land in both WPA and PC sites. The solid line is the g(r) function curve generated based on the ALS segmentation-based tree location point map, and the dashed lines are the 95 % confidential interval (i.e., envelope) of the Monte-Carlo simulation of the CSR distribution (insets a, b, c, and d). In the insets (e) and (f), the generated g(r) function curves for WPA (inset e) and PC (inset f) sites were compared with the 95 % confidential interval (i.e., envelope) produced based on the Monte-Carlo simulation of the soft-core point process model parameterized using the ALS-based metric.

ALS-Determined Parameters of Tree Spatial Distribution Models
Once the tree spatial distribution pattern is determined, an appropriate empirical statistical model can be used to approximate the tree spatial distribution pattern. Based on the coordinates of trees and crown sizes, we obtained the unknown parameters of different statistical models for six natural forest plots in both WPA and PC sites (Table 3). For forest plot WPA-HC with clustered distribution pattern, three unknown parameters including N, rN, and mN were obtained for the Neyman-A model using ALS-(visual-) based data as 0.016 trees / m 2 (0.018 trees / m 2 ), 5.09 m (3.25 m), and 1.37 (1.50), respectively. Moreover, parameters β, HC, rg, γ and s of hybrid Gibbs model were determined as 0.007 (0.009), 1.65 (1.50) m, 3.5 (3.0) m, 1.92 (2.11), and 3 (3) based on the ALS-(visual-) based data, respectively. In terms of the regular distribution WPA-MC and PC-HC, the determined parameters (β, σ, and κ) for the soft-core point process models were 0.02 (0.045), 5.45 (4.50) and 0.50 (0.39) using the ALS-(visual-) based data for plot WPA-MC. In plot PC-HC, the ALS-and visualbased unknown parameters (β, σ, and κ) of the Soft-core point process model were 0.027 and 0.035, 4.45, and 3.70, 0.40 and 0.32, respectively. As for the plots WPA-LM, WPA-LC, and PC-MC were determined as random distribution, the ALS-and visual-based density (λ) parameters were the same as the ones of forest plots WPA-LM (0.006 trees / m 2 ) and WPA-LC (0.007 tree / m 2 ). Moreover, the ALS-and visual-based densities were 0.008 trees / m 2 and 0.009 trees / m 2 for forest plot PC-MC. Table 3. ALS-determined input parameters for the Poisson, soft-core, Neyman-A, and hybrid Gibbs models.

Visual-Based Values
Hybrid Gibbs (WPA-HC) According to the methods described in section 2.5, we determined the unknown parameters of the Soft-core point process models for the complete coverage in both WPA and PC sites based on the ALS-based metrics. The detected number of trees was 2,745 and 7,669 for WPA and PC sites, with the average crown radii as 5.49 m and 5.03 m, respectively. The determined unknown parameters including σ, κ, and β were 4.73 (6.52), 0.60 (0.75), and 0.10 (0.19) for WPA (PC) site (Table 4). Furthermore, for the whole study area of both WPA and PC sites, the ALS-based g(r) function completely fell into the envelope lines generated by the ALS determined soft-core point process model (Figure 8e, f). In addition, we also determined the three parameters (i.e., β, κ, and σ) of modeled soft-core point process in the R programming environment as 0.013, 0.1, and 4.64 for forest Plot-1 and 0.025, 0.1 and 4.50 for forest Plot-2, respectively. The determined parameter of the Poisson point process for random distribution were 0.007 and 0.012 for modeled forest Plot-3 and Plot-4, respectively. In terms of the simulated clustered forest distribution pattern, the three unknown parameters (i. e., N, rN, and mN) were 0.002, 7.578, and 3.083 for forest Plot-5 and 0.002, 4.378, and 3.783 for Plot-6, respectively using the Neyman point process.

Validation of ALS-Determined Tree Spatial Distribution Models
To verify the applicability of the statistical models with determined parameters derived from ALS data we compared the visual-based g(r) function curves with the 95 % envelope lines generated based on the ALS determined statistical models for six natural forest plots in both WPA and PC sites (Figure 9). It was found that a partial visual-based g(r) function curve was out of the envelope lines for forest plot WPA-HC (Figure 9a). The peak value of the ALS-based g(r) function shifted toward the right side, indicating that the tree spatial distribution pattern characterized by ALS data has less clustered strength. The corresponding larger r value showed that the clustering occurred at a larger spatial scale. This discrepancy might be attributed to the errors of individual tree segmentation result from the small tree crowns adjacent to the big tree crowns. The mis-segmentation will result in bigger ALS-based crown radius values. However, the visual-based g(r) function curves fall entirely into the envelope lines generated based on the ALS determined soft-core point process and Poisson model in forest plots WPA-MC (Figure 9b (Figure 9f). It indicated that appropriate statistical models with determined parameters derived from ALS data had been successfully applied to approximate the tree spatial distribution patterns.
As for the forest WPA-HC, we plotted the ALS-and visual-based g(r) function curves with the envelope lines generated based on the ALS determined Neyman-A model (Figure 9g). It was found that both ALS-and visual-based g(r) function curves did not fall into the envelope lines suggesting that the Neyman-A model could not approximate the tree spatial distribution. However, the envelope lines obtained based on the hybrid Gibbs model have similar variations with one produced by the visual-based measurements (Figure 9a). Then compared with the Neyman-A model, the hybrid Gibbs model was more effective at modeling the clustered tree spatial distribution.

Spacing Threshold Effects
By conducting a sensitivity analysis of spacing thresholds on the normalized g(r)norm functions, it was shown that the spacing threshold had obvious effects on ALS-based individual tree segmentation results (Figure 10 a -f and Appendix A. Table A1). As the values of spacing thresholds increasing from 1 m to 5 m for six natural forest plots, the number of ALS-based segmented tree crowns decreased, and average tree crown radii increased regardless of the forest plots densities and types (Appendix A. Table A1). Moreover, with the changing spacing threshold, the predicted tree spatial distribution patterns were not constant with various numbers of the segmented tree. For example, in the plot WPA-HC, the tree number and average crown radii obtained using the visualbased measurements were 159 and 3.72 m, respectively. However, we obtained 252 trees and the crown sizes ranging from 1.03 m to 6.79 m, with an average value of 2.96 m, based on ALS data when the spacing threshold was 1 m. Moreover, the predicted tree spatial distribution patterns were mixed when the spacing thresholds were 1 m or 2 m, which was consistent with the pattern produced by the visual-based measurements. However, the tree spatial distribution pattern in WPA-HC became random when the spacing thresholds were as 3 m, 4 m, or 5 m, respectively (Figure 10a). This variation might be explained by the fact that multiple aggregated trees have been identified as one bigger tree when the spacing threshold increased. Then, the selection of the spacing threshold is a key step in the ALSbased method for characterizing the tree spatial distribution patterns. By comparing the numbers of ALS-and visual-based tree crowns in all forest plots, it was found that the minimum tree crown radius was a good reference to set the appropriate spacing threshold. For example, the visual-based minimum tree crown radius (i.e., 1.62 m) in forest plot WPA-HC was similar to the optimum spacing threshold (i.e., 2 m). However, it is not efficient for the mixed forest plot because the appropriate spacing threshold is different for various crown shapes (Figure 11a, d) when applying the ALS-based tree segmentation algorithm. As shown in Figure 11, the broadleaf was segmented into four individual tree crowns, and the two coniferous crowns were delineated as one tree crown when the spacing threshold was set as 2 m (Figure 11b, e). When we set the spacing threshold as 3 m and 1 m, the broadleaf (Figure 11c) and coniferous (Figure 11f) crowns have been identified correctly, respectively. Considering the complexity of forest, we recommend the use of iterative tree crown segmentation method described in section 2.3 to obtain tree locations and crown sizes information from ALS data, especially in landscape or region scales. Figure 11. Graphs showing the effects of spacing thresholds on the number of segmented tree crowns with different colors for broadleaf (insets a, b, and c) and coniferous trees (insets d, e, and f). Insets (a) and (d) are the raw ALS data for broadleaf and coniferous trees. Insets (b) and (e) are the segmented tree crown results when the spacing thresholds were set as 2 m, and insets (c) and (f) are the segmented tree crown results when the spacing thresholds were set as 3 m and 1 m, respectively. The solid white circles represent the tree center locations, different colors in insets b, c, e, and f represent different segmented tree crowns.

Search Radius Effects
The search radius has limited effects on the final accuracy of tree crown segmentation. By plotting the variations of normalized g(r)norm functions with the search radii changing from 1 m to 5 m with the step of 1 m for six natural forest plots in both WPA and PC sites, it was found that similar variation patterns could be observed for all forest plots in both WPA and PC sites (Figure 10g-l). With the search radii increased from 1 m to 5 m, the number of segmented trees decreased, and average crown radii increased in general. For forest plots in the WPA site, as the search radius increasing, the predicted tree spatial distribution patterns changed accordingly since the search radius affected the result of tree crown segmentation. For example, in forest plot WPA-HC, the number of segmented trees was 189, 137, 129, 113, and 108 as the search radii increasing from 1m to 5 m with the step of 1 m. Furthermore, the average value of crown radii was 3.13 m, 3.38 m, 4.65 m, 6.1 m, and 7.21 m with an increasing search radius (Appendix A. Table A2). Moreover, when the search radii were beyond 4 m, the tree spatial distribution pattern shifted from clustered to random distribution (Figure 10g).
The natural forest plots showed constant distribution patterns such as regular and random forest plots PC-HC and PC-MC (Appendix A. Table A2). Furthermore, through comparing the ALS-and visual-based tree crown radius, the ALS-based tree crown sizes showed better performance in natural forest stands with small averaged standard deviation (i.e., 0.95) of crown sizes compared with the one (i.e., 1.38) in a well-managed urban forest. A thumb rule to determine the search radius was that the search radius should be less than the minimum crown diameter in the plot.

Effect of Overlap on Individual Tree Segmentation
The complexity of forest structure will affect the accuracy of tree crown segmentation. Forest stands with monodominant and separated tree crowns will achieve better accuracy and efficiency compared with the mixed forest with substantial overlapping between tree crowns. The overlapping effect makes it difficult to identify the boundaries of individual tree crowns accurately, especially in adjacent regions. As shown in Table 5, the increasing overlapping distance between tree crowns affected the accuracy of tree crown segmentation results. In the case of two coniferous trees with the tree height and crown diameter of 24.50 m and 8.03 m, it showed that the number of classified points for Tree-1 changed from 2968 to 3757 by changing the overlapped distances between two tree crowns from 0.5 m to 4.0 m. However, the number of points for Tree-2 decreased from 2894 to 2105. Furthermore, the average error of segmented tree crown diameter changed from 0.22 m to 2.12 m as the overlapped distance increasing. It showed that the overlapping distance of tree crowns should be less than 2.50 m to achieve better accuracy (> 80 %) of tree crown segmentation in the coniferous forest in our work. For two broad-leaved trees with tree height and crown diameter of 10.52 m and 10.26 m, the average error of segmented tree crowns increased from 0.42 m to 2.92 m as the overlapping distance increasing from 1 m to 4 m. For the overlapped crowns mixed with conifer (i.e., tree height is 15.38 m and crown diameter is 7.59 m) and broadleaf (i.e., tree height is 10.52 m, and crown diameters is 10.26 m), the average errors of tree crown segmentation increased from 0.44 m to 6.95 m with the overlapping distance increased from 1 m to 3.5 m. It showed that the overlapping range of tree crowns should be less than 2 m to achieve higher accuracy (80 %) of tree crown segmentation in the mixed forest in our work. However, it should be noted that the reasonable overlap rate varied with tree species, crown shape, crown size, relative tree height difference, and canopy cover. In addition, the comprehensiveness of the acquired forest point cloud data greatly affects the accuracy of tree crown delineation. For example, the shadow effects resulted from the tall tree might prevent the penetration of laser beams from capturing the crowns of short trees. This shadow effect could be improved by the cross-direction method or increasing the density of the flight path.
In addition, we compared our segmentation result with the ones obtained from the method proposed by Silva et al. (2016) using liDR package in the R programming environment [58,59]. It was shown that the liDR-based method tended to over-segment.

Effect of the Hard-Core Radius on Tree Spatial Distribution Pattern
The hard-core radii of individual tree crown within a certain range, which vary with tree species, crown sizes, and environmental conditions, usually has relatively stable effects on the tree spatial distribution pattern. Taking the forest plot WPA-HC for example, we compared the visual-based g(r) function curves with the envelope lines of the hybrid Gibbs models parameterized using visual-based hard-core radii in forest plot WPA-HC ( Figure 12). It was found that, overall, the visual-based g(r) function curves could well fall into the envelope lines generated based on the hybrid Gibbs model with the hard-core radii of 1.4 m, 1.5 m, 1.9 m, and 2.0 m. More specifically, when the hard-core radius was 1.4 m, the visual-based g(r) function curve was out of the range of envelope lines when the g(r) < 1 at the regular distribution scale and shifted toward the right of the envelope lines (Figure 12a). It might be attributed to the fact that the smaller hard-core radius (e.g., 1.4 m) used in the hybrid Gibbs model could not explain the true competitive effects denoted by a visual-based hard-core radius (i.e., 1.5 m) well. When the hard-core radii ranged from 1.5 m to 1.9 m, the g(r) function fell within the 95 % envelope lines (Figure12b, c). However, as shown in Figure 12d, the visual-based g(r) function curve was out of the envelope lines produced based on the hybrid Gibbs model parameterized with a larger hard-core radius (i.e., 2 m). It indicated that the hybrid Gibbs model with determined parameters might not be able to capture the aggregation effect among tree crowns in this case. Therefore, the hard-core radius for a given forest plot is not a fixed value, especially in cases with various tree species and crown sizes. The fact that different tree species and crown sizes may have different hard-core radii might explain the unstable radii of hard-core of a forest stand.

Spatial Variations of Tree Spatial Distribution Patterns
The tree spatial distribution patterns usually exhibited varying patterns with the changes in the spatial scales of the observation windows, especially in the forest with a clustered distribution. For example, the random distribution pattern could be found within the observation windows with a smaller size in the natural forest PC site ( Figure 13b). As the size of the observation window increased to 135 m × 135 m, the spatial distribution pattern gradually shifted into regular patterns due to the competitive effects (Figure 13c). The distribution pattern became a clustered one in 500 m × 500 m observation window (Figure 13d). The variations of the distribution patterns with the changing of observation window sizes have been shown clearly in its g(r) function curves. Therefore, an appropriate observation window should be determined before exploring tree spatial distribution patterns according to the interest of the study. A larger observation window should be used if the primary interest is to investigate the effect of environmental factors on tree spatial distribution patterns such as soil fertility. The observation window should be big enough to encompass the typical variations in environmental conditions. Wiegand et al. [17] found similar results and recommended that a suggested number of tree crowns was more than 70 in determining the size of observation windows. In this study, a smaller observation window was preferable to minimize the effects of environmental factors and focus on the potential interactions between tree crowns. Moreover, the tree spatial distribution patterns are not constant but vary with human practices such as deforestation or afforestation, disturbances such as fire or insect damages, or natural regeneration growth. It is a direction deserving of further research.

Tree spatial Pattern Effects on the Probability of Viewing Ground Area
Tree spatial distribution patterns showed obvious effects on the probabilities of viewing the ground surface area from a range of viewing directions. By following an existing method [55], we first obtained the probabilities of viewing ground surface area in five directions for six artificially built ALS forest plots. For three forest plots with the same number of trees, the probabilities of viewing ground surface area increased from 55.63 %, 59.52 % to 66.01 % as the tree spatial distribution patterns changed from regular (Plot-1), random (Plot-3) to clustered (Plot-5) distribution in nadir direction (Table 6). Similar changing trends of viewing ground probabilities were observed for other viewing directions in principle plane.
Moreover, the determination accuracy of unknown parameters for a given tree spatial distribution pattern showed great effects on the probabilities of viewing the ground surface. The same tree spatial distribution pattern might show different viewing probabilities by modifying the parameters of a given statistical model. For example, as shown in Table 6, the values of parameter λ were 0.007 and 0.012 in the randomly distributed forest Plot-3 and Plot-4. Apparent changes were observed from 59.52 % to 39.92 %, from 56.87 % to 36.88 %, from 36.15 % to 16.67 %, from 47.06 % to 26.79 % and from 30.37 % to 12.99 % at five different viewing directions in the principle plane. Similar variations could also be observed for regular (i.e., forest Plot-1 and Plot-2) and clustered (forest Plot-5 and Plot-6) distributions as well. These changes have been shown great effects on the BRDF modeling of a forest canopy [60,61]. It was shown that both the selection of statistical models and the determination of unknown parameters are important in quantitatively characterizing tree spatial distribution patterns. The ALS-based method in this study provides an objective and effective way to quantitatively determine tree spatial distribution patterns compared with the traditional aerial photography-based methods [16,33].

The Effects of Tree Location Proxy Selection
Appropriate proxy selection of tree location is essential for tree spatial distribution pattern studies. Both tree stem and treetop point location could be used to represent individual tree locations. In terms of the tree spatial distribution pattern, there should be no apparent differences between tree spatial distribution patterns if the treetop points coincide with stem locations in a forest stand. While different spatial distribution patterns could be observed using the treetop and stem location-based approaches when treetop points are off from stem locations due to competition between tree crowns for light and various resources [47]. For example, a forest stand with random stem distribution pattern might produce clustered treetop distribution pattern. In this case, the probability of viewing the ground surface from the sky in GO models [11,62] could be better characterized quantitatively by the distribution pattern of tree crowns instead of tree stems due to the occlusion effects and shadows between tree crowns.
In addition, the conventional field-based stem location method is labor-intensive, timeconsuming, and high cost, which limits its application from achieving actual tree spatial distribution patterns at the landscape or regional levels. Then, the motivation of current work was to help the modeler to determine the unknown parameters of statistical methods used in their GO models based on the ALS-derived treetop location and crown size information with high efficiency and accuracy and low-cost. The more accurate tree spatial distribution pattern will further improve the simulated directional forest canopy reflectance and the retrieval accuracy of forest structural parameters such as LAI based on the remotely sensed spectral information.
In summary, we first need to determine the size of the observation window according to our research goal. Second, we can characterize tree spatial distribution patterns using the proposed method in the present study based on the individual tree segmentation results. It should be noted that the segmentation result will affect the result of quantifying tree spatial distribution. Currently, many methods have been proposed to delineate individual trees using ALS data for different forest types [41,[63][64][65][66]. Thus, an appropriate segmentation method can be chosen in terms of forest structure characteristics in the study area. In addition, the multiple iterations of different segmentation methods are recommended at the regional scale, especially with high tree species richness. Furthermore, how to delineate tree crown in the heterogeneous forest with heavily overlapped tree crowns requires further investigation.

Conclusions
In this study, we developed an ALS-based method to characterize quantitatively different types of tree spatial distribution patterns. Based on our results, we concluded that: (1) the proposed ALSbased approach provides a practical and effective way to determine the unknown parameters of various statistical models for characterizing tree spatial distribution patterns purpose; (2) hard-core radius can quantify the competitive effects between tree crowns effectively using the proposed individual tree crown model; (3) the Poisson and soft-core models are capable of describing random and regular distributions, respectively. For clustered distribution, compared with the Neyman-A model, the hybrid Gibbs model can better quantify the tree spatial distribution patterns across various spatial scales. This work will be beneficial to the studies of radiative transfer modeling and spatiotemporal variations of forest canopy BRDF.
Author Contributions: G.Z. and X.W. conceived and designed the experiments; X.W. performed the experiments; X.W. analyzed the data; L.M.M. contributed data; X.W. wrote the paper; Z.Y. contributed analysis tools; G.Z. modified the paper. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.