Point process models for mixed sessile forest stands

Abstract

• Background

Growth modelling of complex stands calls for the use of spatially explicit single-tree models. Such models require spatially explicit tree locations as the initial state to run simulations. Given the cost of such data, virtual forest stands, where tree locations are simulated, are generally used as the initial state.

• Purpose

The purpose of this study was to present models for simulating the spatial structure of complex stands. It focused on mixed oak–Scots pine stands of the Orleans forest (France) and on the spatial structure of canopy trees.

• Methods

The spatial structure of the oak–pine stands was modelled with appropriate point process models. The models consisted of a combination of Poisson, Neyman–Scott and Soft core. Simulation of the point process models was based on precise characterisation of the studied stands. Twenty-five 1-ha oak–pine plots were characterised by the Ripley function. The models were then fitted to the identified spatial structure to reproduce univariate and bivariate spatial patterns in each spatial type.

• Conclusion

This paper provides an approach for general modelling of a spatial structure of a particular mixture and may be enriched by other point process models for other types of mixed stand.

1 Introduction

Complex forest stands, of mixed or uneven ages, are an important silvicultural regime for sustainable forest management. As the area of such stands increases, due to ecological, social and economic demands (Buongiorno et al. 1995; Caquet et al. 2010), mixed forest management is faced with some research issues. Their dynamics are more difficult to elucidate than those of pure or even-aged stands. More particularly, in a context of changing environmental and social conditions, it is urgent to develop relevant stand description tools and models of stand dynamics to predict the evolution of such forests (Dobbertin 2009).

Modelling the dynamics of mixed and uneven-aged stands is thus a major challenge in forest research. Models, particularly growth models, and their simulations are useful not only for defining management patterns and silvicultural guides but also for evaluating current and future resources (Biging and Dobbertin 1995; Pacala et al. 1996).

In growth modelling, a distinction is made between three categories of growth models, based on their spatial resolution: single-tree level, diameter class level and stand level. Of these, the single-tree approach seems most appropriate for understanding stand growth as affected by competition between individuals of different species (Gourlet-Fleury and Houllier 2000; Courbaud et al. 2001; Goreaud et al. 2007). Such models are spatially explicit and simulate the evolution of each individual tree depending on its own characteristics and its local environment.

However, forest researchers and managers cannot easily implement these models because they require spatially explicit tree locations to run simulations. There is therefore an increasing demand for spatial tree data sets for the models (Pommerening 2006). Such data are scarce and costly, and it is often necessary to simulate the necessary tree positions in a virtual stand. Consequently, using virtual stands as the initial state is a usual alternative. It appears important for virtual stands to be realistic, i.e. statistically as similar as possible to real forest patterns.

The aim of this study was to present a model for simulating the spatial structure of mixed/uneven-aged forest stands. A clear understanding and simulation of spatial dependence and its quantification are important aspects of forest research and sustainable management of forests (Goreaud 2000; Pommerening 2006). Indeed, tree growth is a reaction to spatial context and, conversely, growth processes affect spatial structure.

This modelling approach was applied to mixed oak–Scots pine stands of the Orleans forest (France). It was based on a previous study (Ngo Bieng et al. 2006) presenting the various spatial organisations observed in these forest stands. This paper describes the method used in modelling real spatial structures identified in the previous study, along with the results of the modelling approach, and discusses the merits of such models of spatial structure.

2 Methods

2.1 Study area

2.1.1 The stands

The study area was located in the Orleans state forest, one of France’s largest public woodlands (35,000 ha) in central France (47°50′S, 2°30′E). It was dominated by oak and Scots pine, but also included trees of various other species, such as birch (Betula pendula), wild service tree (Sorbus torminalis) and hornbeam (Carpinus sp.). The study focused on the southern part of the forest, where we were specifically studying mixed stands of sessile oaks (Quercus petraea) and Scots pine (Pinus sylvestris).

2.1.2 Sample plots and data

We set up 25 square (100 × 100 m) or rectangular (50 × 100 m) plots (0.5 or 1 ha) with uniform site conditions and containing more than 10% of mixed species canopy. The 25 plots were selected to reflect a range of mixture rates and tree size distributions (unpublished data).

The 25 plots were mapped to study the spatial structure of mixed oak–pine stands. Inside each plot, the exact location of each tree (diameter at breast height >7.5 cm) was measured using a theodolite. For each located tree, we recorded the species and circumference at breast height (at 1.3 m) and whether it was part of the canopy tree subpopulation or the understorey tree subpopulation (see the following paragraph).

2.1.3 Defining subpopulations

The 25 plots showed a wide range of mixture rates, heights and circumferences, prompting us to define various subpopulations corresponding to homogeneous groups of trees (Ngo Bieng et al. 2006). The defined subpopulations revealed different species (represented by at least ten individuals per hectare). The subpopulations were also based on vertical strata corresponding to canopy trees and understorey trees. This definition of subpopulations took into account the vertical gradient of light availability and also, implicitly, the size classes (characterised by high variability in tree diameters); the largest in terms of diameter were in the canopy.

This study focused on modelling spatial interactions between trees in the canopy (with an estimated height of 18 m or more). We considered two subpopulations in this work: oak canopy trees and pine canopy trees as the other species were not significantly represented in the canopy of each plot.

2.2 Spatial structure analyses

Ripley’s K function (Ripley 1977; Diggle 1983; Cressie 1993) was used to study the univariate pattern, i.e. the spatial pattern of each subpopulation. This function has been used to examine spatial point patterns in forest ecosystems (Perry et al. 2006) and quantifies the degree of clumping and/or overdispersion of the pattern on multiple scales. Ripley’s K function is a function of the mean number of other trees found within a distance r from a typical tree and defines different degrees of random, clustered or regular spatial organisation.

The spatial interactions between two subpopulations, i.e. the bivariate pattern, were investigated by the K ₁₂(r) intertype function. This function is based on the same principle as the K(r) function, but is used to characterise the relationship between two subpopulations (1) and (2) in the same plot. It defines different degrees of interspecific independence, attraction or repulsion (Lotwick and Silvermann 1982).

As recommended by Besag (in Ripley 1977), we used L(r) and L ₁₂(r), square-root transformations of K(r) and K ₁₂(r) functions (Ngo Bieng et al. 2006).

2.3 Point process models

When analysing the spatial structure of a forest stand, the stand is considered to be a finite set of trees T _i (i = 1,..., n) in a bounded area, characterised by their locations X _i  = (x _i , y _i ) and their marks m _i (Stoyan and Penttinen 2000). It is thus appropriate to apply a planar (marked) point process method: the tree locations are “points” and the tree characteristics (species, diameter at breast height or dbh, etc.) are “marks”.

A spatial point process is a stochastic model that governs the location of points in space (Diggle 1983; Cressie 1993; Illian et al. 2008). The relevance of point processes for the study and simulation of forest stand spatial structures is largely documented. Tomppo (1986), Stoyan and Penttinen (2000) and Comas and Mateu (2007) reviewed the main point process methods applied in forestry. In forestry applications, we are interested in the spatial point pattern of trees in a forest stand as a result of a given point process. Many point processes have been used to generate stand spatial structures. In this paper, the point processes used were:

The homogeneous Poisson process, which simulates a random pattern, i.e. the configuration of independently and uniformly distributed points. It is a null or reference model and is the ideal model for complete spatial randomness (CSR; Tomppo 1986; Comas and Mateu 2007).

The Neyman–Scott (NS) process, which is a Poisson cluster process, was designed to generate clustered or clustered spatial point patterns (Tomppo 1986).

Simple sequential inhibition is a “hard core” process which simulates a univariate pattern of regularity and a bivariate pattern of repulsion (Cressie 1993; Tomppo 1986). It is a pairwise interaction process where there are no pairs of points that are closer than a threshold distance or “hard core” distance. In our modelling approach, we used “soft core” (SC), which is based on the same principle as hard core, but with an interpoint interaction decreasing smoothly with distance.

These point processes were combined in the simulation procedure in order to reproduce accurately the average spatial arrangement of each subpopulation within each identified spatial type of the studied stand: one point process for each subpopulation, and the two processes selected enabled realistic simulation of the spatial structure for each subpopulation, while taking into account the intertype structure between the two subpopulations.

First of all, we tested a Poisson forest to assess whether the point configurations in a spatial type could be regarded as generated by a Poisson process.

2.3.1 M₀: the null model consisting of a combination of two Poisson processes

This model consisted of the simulation of two point patterns (one for each species), with a random spatial structure and an intertype structure of independence. The two point patterns, corresponding to the number of oaks n ₁ and the number of pines n ₂, were simulated by two conditional Poisson processes (conditioning on n ₁ and n ₂, respectively) over a given area. This restriction of the Poisson process to the bounded area, conditional a number n, yields to a new point process which is actually the binomial point process (Illian et al. 2008). The tree positions were randomly distributed in each point pattern. This model M₀ was the representation of the null hypothesis of random structure for each specific spatial structure and independence for the intertype structure.

We then tried more complicated combinations in accordance with identified spatial characteristics.

2.3.2 M₁: a combination of the Neyman–Scott and soft core processes

In order to reproduce interactions between trees (intraspecific regularity and interspecific repulsion), simulation was carried out on the basis of sequential inhibition processes: trees were placed sequentially in a bounded region, and a new position was accepted if its distance from all previously accepted trees was at least a threshold distance of regularity or repulsion.

Consequently, the M₁ model consisted of the simulation of two point patterns: the first point pattern, corresponding to pine, was simulated with a combination of the Neyman–Scott and SC processes. The pine points were simulated in clusters by the NS process, with an additional constraint of regularity between simulated points provided by an SC process. Trees were placed sequentially within clusters, and a new position was accepted if its distance from all previous pines was less than a given threshold distance of regularity (d _reg). After simulation of the pines, the oak points were simulated by an SC process. Oaks were placed sequentially in the simulated area, and a new tree was accepted if its distance from all the pine points was at least a given threshold distance of repulsion (d _rep).

With the SC process, interactions are softer than the interaction defined in a simple sequential inhibition process. Indeed, the interaction patterns of trees (regularity and repulsion) are often softer in real forest stands. To simulate such a soft interaction process, during simulation, we accepted the positions of trees simulated at a distance below the given threshold distance, with a probability depending on the distance. This probability p (simulated distance/threshold distance) varied linearly from p = 0 at a null distance to p = 1 at the given threshold distance. Consequently, the regularity was not hard, but soft: there may have been some pine points at a distance of less than d _reg from the other pine points. Similarly, repulsion was not hard: there may have been some oak points at a distance less than d _rep from the pine points.

Lastly, the M₁ model had four parameters, given n ₁ oaks and n ₂ pines in a given area: the number of clusters of pines (n _ag), the radius of the clusters (r _ag), the minimum distance between pines (regularity distance, d _reg), and the minimum distance between pines and oaks (repulsion distance, d _rep).

2.3.3 M₂, a variation of the M₁ model combining NS and SC processes

Model M₂ was similar to the previous model. The pine points were first simulated with a Neyman–Scott point process, with an additional constraint of soft regularity between the points. The oak points were then simulated in soft repulsion with the pine points. The difference from model M₁ was that, here, the probability p of accepting an oak point closer to a pine than the threshold repulsion distance d _rep was constant and did not vary with the distance.

This model thus had five parameters: the same four parameters as the previous model, plus a fifth parameter, p, the constant probability of a pair intertype configuration being accepted regardless of d _rep.

2.4 Fitting procedure

The models presented above were used to simulate a bivariate point pattern in a defined area, a point pattern for oak and a point pattern for pine. Once the simulated point pattern was obtained, its spatial structure could be described with the Ripley and intertype functions: the spatial structure of the simulated oak point (L ^sim_o (r)), the spatial structure of the simulated pine point (L ^sim_p (r)), and the intertype structure between the simulated oak point and the simulated pine point (L ^sim_op (r)).

The parameters of these models could then be fitted so that the Ripley and intertype simulated curves were statistically as similar as possible to the Ripley and intertype mean curves of an identified spatial type (L ^real_o (r), (L ^real_p (r), (L ^real_op (r)). Following the minimum contrast method (Diggle and Gratton 1984), we used the least squares criterion (LSC) defined in Eqs. 1 and 2 to minimise the distance between the simulated and real curves. We wanted to reproduce not only the mean spatial characteristics in each identified spatial type but also the individual variability inside each type. Thus, the least squares criterion was computed on two levels:

• (1) the squared Euclidean distance between the real and simulated mean curves (Eq. 1);

  $$ {\hbox{LS}}{{\hbox{C}}_1} = \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\rm{mean}}^{\rm{sim}}}_{{{\rm{(o)}}}}}(r) - L{{_{\rm{mean}}^{\rm{real}}}_{{{\rm{(o)}}}}}(r))}}\nolimits^2 } \right]} + \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\rm{mean}}^{\rm{sim}}}_{{{\rm{(p)}}}}}(r) - L{{_{\rm{mean}}^{\rm{real}}}_{{{\rm{(p)}}}}}(r))}}\nolimits^2 } \right]} + \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\rm{mean}}^{\rm{sim}}}_{{{\rm{(op)}}}}}(r) - L{{_{\rm{mean}}^{\rm{real}}}_{{{\rm{(op)}}}}}(r))}}\nolimits^2 } \right]} $$

  (1)

• (2) the squared Euclidean distance between the real and simulated standard deviations (Eq. 2).

  $$ {\hbox{LS}}{{\hbox{C}}_2} = \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\sigma }^{\rm{sim}}}_{{{\rm{(o)}}}}}(r) - L{{_{\sigma }^{\rm{real}}}_{{{\rm{(o)}}}}}(r))}}\nolimits^2 } \right]} + \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\sigma }^{\rm{sim}}}_{{{\rm{(p)}}}}}(r) - L{{_{\sigma }^{\rm{real}}}_{{{\rm{(p)}}}}}(r))}}\nolimits^2 } \right]} + \sum\limits_{{r = 2}}^{{30}} {\left[ {\mathop{{(L{{_{\sigma }^{\rm{sim}}}_{{{\rm{(op)}}}}}(r) - L{{_{\sigma }^{\rm{real}}}_{{{\rm{(op)}}}}}(r))}}\nolimits^2 } \right]} $$

  (2)

The real mean curves and standard deviations were obtained from the values in each spatial type of the Ripley and intertype functions of the plots that formed a spatial type. The simulated mean curves and standard deviation were obtained from the values of 1,000 runs of a model with a given set of parameters.

We were interested in the parameters that minimised the difference between the model and the empirical data, and so, a selection was made from different parameterisations. The least squares criteria were minimized by searching over the parameter space of each model. The parameter space was as large as possible, for exampl,e for the M₂ model: n _ag from 25 to 45, r _ag from 5 to 15, d _reg from 1 to 8, d _rep from 5 to 13 (every n + 1 for the fourth parameter), and p from 0 to 0.5 (every n + 0.05 for this parameter).

The simulation and fitting procedures were implemented in JAVA, and we built a programme for the analysis of spatial structure dedicated to analyses and simulations of forest spatial structures (Ngo Bieng 2007). The point process models were also implemented and available on the CAPSIS simulation platform (De Coligny et al. 2004).

3 Results

3.1 Summary of the spatial types identified in the studied stands

The canopy spatial structure of the studied stands, identified with the L(r) and L ₁₂(r) functions, highlighted variability in the spatial pattern of the plots: for the univariate pattern, a gradient ranging from significant aggregation on different scales to randomness of the main two species; for the bivariate pattern, the intertype structure also ranged from clear interspecific repulsion to independence.

The identified spatial arrangements organised in a typology of spatial structure were presented in Ngo Bieng et al. (2006). A spatial type consisted of plots with similar spatial structure. Fig. 1 shows the spatial characteristics of one of the identified spatial types corresponding to the spatial structure of young stands whose dynamics and evolution are of interest for managers. This type was characterised by a structure not differing significantly from randomness for oaks (Fig. 1a) and slightly clustered for pines (Fig. 1b), with a tendency to regularity of pines at small distances. The intertype structure (Fig. 1c) ranged from interspecific independence to interspecific repulsion at small distances. In what follows, the point process models were fitted to reproduce the spatial structure of that type.

Fig. 1

Modified L(r) functions for oaks and pines in the canopy for the spatial type presented: a oaks Lo(r), b Scots pines Lp(r) and c intertype Lop(r) functions. The thin lines with different symbols identify the spatial characteristics of each of the six plots. The thick lines are the mean curves of this type

3.2 Results of the fitting procedure

3.2.1 M₀: the null model

We fitted this model to the data by simulating two conditional Poisson processes: the first process conditioned on n ₁ oaks per hectare, and the second process conditioned on n ₂ pines per hectare; n ₁ (119 oaks) and n ₂ (123 pines) were obtained from real data (average number of oaks and pines over the six plots with the spatial type presented). We simulated 100 runs of this null model to compare them with the empirical data. We obtained a squared Euclidean distance of 49.6 (LSC₁ = 34.7 between the real and simulated mean curves and LSC₂ = 14.9 between the real and simulated standard deviations).

Fig. 2 presents a graphic comparison of the mean L(r) and intertype functions, empirical in black and simulated in grey. A simulated 95% confidence envelope was built from 100 runs of the null model (pale grey lines).

Fig. 2

Empirical (black line) and simulated by the null model (grey line) L(r) and L ₁₂(r) functions, left to right a L _o(r), b L _p(r), c L _op(r); dashed lines indicate 95% simulated confidence envelopes

3.2.2 M1: a combination of the Neyman–Scott and soft core processes

This model consisted of the simulation of n ₂ pines with a combination of the NS and SC processes (to monitor regularity) and the simulation of n ₁ oaks with an intertype SC process. This model was fitted to our data using the fitting procedure presented above. The fitted parameters for model M₁ were: n _ag = 23, r _ag = 8, d _reg = 10 and d _rep = 4, corresponding to a total least squares criterion of 19.47 (LSC₁ = 10.07; LSC₂ = 9.4).

Fig. 3 presents a graphic comparison of the mean L(r) and intertype functions, empirical in black and simulated in grey. A simulated 95% confidence envelope was built from 100 runs of the model M₁ (pale grey lines).

Fig. 3

Empirical (black line) and simulated by the model M₁ (grey line) L(r) and L ₁₂(r) functions, left to right a L _o(r), b L _p(r), c L _op(r); dashed lines indicate 95% simulated confidence envelopes

3.2.3 M₂, a variation of the M₁ model combining NS and SC processes

Like model M₁, this model consisted of the simulation of n ₂ pines with a combination of the NS and SC processes and the simulation of n ₁ oaks with an intertype SC process. Here, there was one more parameter, p, which was the constant probability of a pair intertype configuration being accepted regardless of d _rep. This model was fitted to our data in the same way as previously, and the fitted parameters for the model were: n _ag = 38, r _ag = 8, d _reg = 10, d _rep = 4 and (v) p = 0.15, corresponding to a total least squares criterion of 17.7 (LSC₁ = 5.8; LSC₂ = 11.9).

Fig. 4 presents a comparison of the mean L(r) and intertype functions, empirical in black and simulated in grey. A simulated 95% confidence envelope was built from 100 runs of the model M₂ (pale grey lines).

Fig. 4

Empirical (black line) and simulated by the model M₂ (grey line) L(r) and L ₁₂(r) functions, left to right a L _o(r), b L _p(r), c L _op(r); dashed lines indicate 95% simulated confidence envelopes

4 Discussion

With an LSC of 49.6, model M₀ seemed inappropriate for simulating the spatial type presented. Model M₀ reproduced the specific spatial structure of oak quite satisfactorily: the mean empirical L(r) curve fell within the simulated confidence envelope. The results were worse for the specific spatial structure of pine and the intertype structure. Not surprisingly, the simulated pattern of pine was more random, and the empirical L(r) was above the simulated confidence envelope. The empirical L ₁₂(r) fell below the confidence envelope. Model M₀ did not reproduce the identified interspecific repulsion.

Models M₁ (LSC = 19.47) and M₂ (LSC = 17.7), where there were interactions between points, reproduced the presented spatial structure more successfully than the null model M_0. The empirical L(r) and L ₁₂(r) curves fell within the simulated confidence envelopes. In particular, a comparison of the empirical and simulated mean curves showed that the simulations with the M₂ model were closest to the empirical curves. The M₁ model simulated a more clustered spatial structure of oak and a more random structure of pine compared with the real curves. For both models, the simulation minimised repulsion on small scales.

The point process models for mixed oak–pine stands, consisting of a combination of appropriate point processes, were based on a typology of spatial structure. The point process models therefore reproduced the spatial structure identified within each type of the typology. This modelling approach, illustrated for one identified type, was reproducible for all the other canopy spatial types. However, it implied different combinations of point processes (Ngo Bieng 2007). For example, as highlighted in Section 3.1, the typology of spatial structure defined different spatial patterns for the plots, consisting in a gradient from mono-specific clusters in interspecific repulsion to a more “mixed” interspecific pattern. The “clusters in repulsion” pattern corresponded to older stands, with thinned clusters due to thinning from diameter and crown growth, competition for resources or sometimes management (Moeur 1997; Batista and Maguire 1998). This regular spacing in adult clusters resulted in the simulation of a “hard” regularity with the hard core process, and the more clustered pattern resulted in different fitted parameters of a double Neyman–Scott process (for oak and pines, smaller number (n _ag) of larger clusters (r _ag)).

The least squares criterion was used to fit the point process models to the data. The set of parameters was chosen in a large space, from which the models minimised the squared Euclidean distance between the empirical and simulated curves. This approach corresponded to a traditional method of parameter estimates in point process statistics: the minimum contrast. It has been used by some authors (Batista and Maguire 1998; Kokkila et al. 2002), who performed the fitting procedure by a comparison with real and simulated curves. In our case, we fitted the models by comparing real and simulated mean curves. This fitting procedure defined parameters describing salient characteristics of the stands (number of clusters per hectare, number of trees per cluster, etc.). The real and simulated standard deviations were also taken into account when computing the squared Euclidean distance. Taking into account the standard deviations when fitting the models was original, and it enabled us to allow for the real variability within each spatial type. However, this individual variability was only obtained with a small number of plots.

The fitting results were satisfactory, the empirical L and L ₁₂ functions could be regarded as a realisation of models M₁ or M₂. The comparison was not only on visual aspects, the simulated envelopes from the fitted models included the L and L ₁₂ functions, but also on the squared Euclidean distance between the empirical and the simulated curves. With the LSC obtained with the M₁ and M₂ models, we were able to conclude that (1) the M₁ and M₂ models gave better results than M₀, and (2) this comparison enabled a quantification of the difference between simulations and data and gave an objective representation of the fitting procedure.

Future prospects may involve:

• Working on the difficulty in accurately reproducing spatial characteristics on small scales. This difficulty has been highlighted in many studies when simulating real stand spatial structures (Moeur 1997; Batista and Maguire 1998; Mateu et al. 1998; Pélissier 1998; Parrott and Lange 2004; Neeff et al. 2005). In our study, this problem especially arose when simulating the intertype structures of repulsion, where there was a minimisation of repulsion on small scales. It may be interesting to give small scales greater weight when fitting models.

• Modelling the spatial pattern of understorey trees and their spatial interaction with canopy trees. Furthermore, it will be interesting to consider quantitative marks (as DBH) and to model marked point patterns for such stands.

• Evaluating the fitted point process models. This evaluation could be made by comparing the simulated and real spatial structures with other spatial summary functions or indices (Stoyan and Penttinen 2000). The evaluation could also be made with data independent of those used for the fitting procedure to see whether point process models can reproduce the spatial structure of other mixed oak–pine stands.

Finally, the accuracy of stand spatial structure simulations can be checked by the response of the spatially explicit single-tree growth model used. It will be interesting to look at the consequences of varying initial spatial structures in growth simulations. Various initial states can be tested, even the CSR model, to assess the gain in precision when using spatial structure models. Errors in spatial pattern can be accepted if they have little effect on the overall growth characteristics of stands (Kokkila et al. 2002).