全 文 ：第 49 卷 第 5 期
2 0 1 3 年 5 月
林 业 科 学
SCIENTIA SILVAE SINICAE
Vol. 49，No. 5
May，2 0 1 3
doi:10．11707 / j．1001-7488．20130515
Received date: 2012 － 08 － 06; Revised date: 2012 － 10 － 12．
Foundation project: Supported by the Scientific Research Funds for Forestry Public Welfare of China(201004026) ; Ministry of Education“Overseas
Experts and Scholars”project．
* Corresponding author: Li Fengri．
树木位置空间模式建模与预测*
金星姬1 李凤日1 贾炜玮1 张连军2
(1. 东北林业大学林学院 哈尔滨 150040; 2. 美国纽约州立大学环境科学和林业学院 锡拉丘兹 NY13210)
摘 要: 传统的林分生长与收获模型不能为林分和生态系统管理提供准确的空间信息。采用 3 种配对势函数的
Gibbs 点过程模型对美国东北部 50 块云冷杉样地里树木的空间分布模式进行模拟。该 Gibbs 模型能够较好地模拟
这 50 块样地中的 82% ～ 84%，但其对完全随机分布和规则分布的样地模拟比对聚集分布的样地模拟效果要好。
使用常用的林分变量如林分密度、公顷断面积、林分平均胸径、平均树高、平均冠幅和冠长建立经验回归模型对
Gibbs 模型的 2 个参数进行预测。结果表明这些回归模型对 81%的样地可以得到满意的模拟效果，其中，100%的
完全随机分布样地、71%的规则分布样地和 56%的聚集分布样地模拟效果较好。选择 3 块样地对树木的模拟空间
位置和实际观测位置的相似性进行对比和说明。
关键词: 空间点模式; 空间点过程; Gibbs 模型; Ripleys K-function; 马尔可夫链 Monte Carlo ( MCMC );
Metropolis 算法; Takacs-Fiksel 方法
中图分类号: S711 文献标识码: A 文章编号: 1001 － 7488(2013)05 － 0110 － 11
Modeling and Predicting Spatial Patterns of Trees
Jin Xingji1 Li Fengri1 Jia Weiwei1 Zhang Lianjun2
(1. School of Forestry，Northeast Forestry University Harbin 150040;
2. College of Environmental Science and Forestry，State University of New York ( SUNY-ESF) Syracuse，NY13210，USA)
Abstract: Traditional forest growth and yield models have been criticized for their inability to provide precise spatial
information for forest and ecosystem management． In this study the spatial patterns of trees within 50 spruce-fir plots in the
Northeast，USA were modeled by a Gibbs point process model with three pair potential functions． In general，82% － 84%
of these 50 plots were modeled well by the Gibbs model． However，the complete spatial random (CSR) and regular spatial
patterns were modeled better than the clustered plots． Further，empirical regression models were developed to predict the
two parameters of the Gibbs model using the available stand variables as predictors such as stand density，basal area，
mean tree diameter，mean tree height，mean crown length，and mean crown width． The simulation results showed that
81% of the 50 plots were satisfactorily predicted by the empirical regression models． Among them，100% of the CSR
plots，71% of the regular plots，and 56% of the clustered plots were predicted well by the empirical regression models．
Three example plots were selected to illustrate the similarity between simulated tree locations and observed ones．
Key words: spatial point pattern; spatial point process; Gibbs model; Ripleys K-function; Markov Chain Monte Carlo
(MCMC); Metropolis algorithm; Takacs-Fiksel method
1 Introduction
Spatial patterns of tree locations in a forest stand
reflect the complex historical and environmental mosaic
imposed by initial establishment patterns，
microenvironment differences， climatic factors，
competing vegetation，management activities，and the
chance success of different species or individuals over
time． The spatial distributions of trees strongly affect
competition between neighboring trees，size variability
and distribution， growth and mortality， and crown
structure ( Kuuluvainen et al．，1987; Kenkel，1988;
Miller et al．， 1989; Hara， 1992; Moeur， 1993;
Rouvinen et al．， 1997; Dovciak et al．， 2001 ) ．
第 5 期 金星姬等: 树木位置空间模式建模与预测
Furthermore， the spatial patterns at tree-scale have
large influences on system-level and community-level
properties of ecosystem functions (Pacala et al．，1995;
Friedman et al．，2001 ) ． In general，positive spatial
dependence among neighboring trees is mainly due to
the effects of the micro-site， whereas inter-tree
competition tends to create negative spatial dependence
among neighboring trees (Fox et al．，2001) ．
Spatial point pattern analysis is commonly used to
study the distributions of the horizontal spacing among
trees within a region of interest． The objective of
spatial point pattern analysis is to measure how
individual trees are located with respect to each other
by using quantitative characteristics and to describe the
“laws” regulating the locations with mathematical
models ( Diggle， 1983; Tomppo， 1986; Pretzsch，
1997 ) ． There are three fundamental spatial point
patterns: complete spatial randomness ( CSR )，
regularity，and clustering． This classification is over-
simplified， but useful at an early stage of spatial
analysis． It is also helpful to formulate a model for the
underlying spatial point process ( Diggle，1983 ) ． A
spatial point pattern is the realization of a spatial point
process ( Cressie， 1993 ) ． Once the spatial point
patterns are recognized，it is desirable to fit appropriate
point process models to the spatial patterns in order to
understand the mechanisms generating the patterns and
to perform analyses of goodness-of-fit for the models．
Detailed review on spatial point process models can be
found in the literature ( Cressie， 1993; Tomppo，
1986; Stoyan et al．，2000) ．
The spatial point process generating CSR is called
a homogeneous Poisson process，such that the number
of events realized by this process 1) follow a Poisson
distribution with mean equal to λ A ，where λ is the
density or intensity of events ( e． g．，trees) and A is
the area of the region A; and 2 ) are independent
random samples from the uniform distribution on A
( Cressie，1993; Diggle，1983 ) ． Intuitively， CSR
means that events are equally likely to occur anywhere
within A and that events do not interact with each other
( Cressie， 1993 ) ． The homogeneous Poisson point
process usually serves as a null model ( Stoyan et al．，
2000) ． If the null model is rejected，more complicated
models are necessary to be tested (Tomppo，1986) ．
One such model is the Gibbs point process model
( or pairwise interaction process model )， which
assumes the interaction between two events ( e． g．，
trees) depends on the relative location of them． It
originates from statistical physics ( Wood， 1968;
Landau et al．，1980) and has been intensively used in
spatial statistics ( Stoyan et al．，1987; Ogata et al．，
1981; 1984; 1985; Diggle et al．，1994) ． The Gibbs
point process model with pair potential functions is a
class of models for point patterns exhibiting interactions
between the events． These models have been applied
in forestry and ecological studies due to their great
flexibility and realism in modeling different spatial
point patterns (Tomppo，1986; Penttinen et al．，1992;
Stoyan et al．，1998; Srkk et al．，1998; Stoyan et al．，
2000) ．
Different methods have been developed to estimate
the parameters of the pair potential functions． Diggle et
al． (1994) summarizes three general methods，namely
approximations of maximum likelihood， maximum
pseudo-likelihood，and Takacs-Fiksel method． These
methods can be routinely used in applications and there
are no restrictions on the parametric form of the pair
potential functions ( Tomppo，1986; Stoyan et al．，
1987; Jensen et al．，1991) ．
Spruce-fir stands in the Northeast， USA are
traditionally important resources for timber production．
Forest managers and researchers have developed a
variety of tools for managing these stands． These tools
include stocking charts and stand density management
diagrams for the pure and mixed softwood stands in the
region ( Solomon et al．，2002 ) ． However，very little
spatial information is available for these stands to date．
The objectives of this research were 1) to explore the
spatial structures of the spruce-fir softwood stands in
the Northeast，USA，2) to fit the Gibbs point process
model with three pair potential functions to the spatial
patterns of trees within the stands，and 3) to develop
empirical regression models for predicting the
parameters of the Gibbs point process model using
available stand variables．
2 Data and methods
2. 1 Data
Fifty ( 50 ) plots were collected from the even-
111
林 业 科 学 49 卷
aged spruce-fir forests in northwestern Maine，USA，
located in the region between 69 #W and 71 #W in
longitude，between 45 #N and 46. 5 #N in latitude，and
between 750 m and 1 200 m in elevation
(Kleinschmidt et al．，1980) ． The plot size ranged from
0. 002 5 to 0. 02 hm2 ． The shapes of the plots were
mostly square with a few rectangular． In each plot，the
position of each tree over 1. 37 m in height was mapped
( Li et al．，2007 ) ． In these plots balsam fir ( Abies
balsamea) and red spruce ( Picea rubens) accounted
for about 95% of total number of trees and 94% of
total volume． Other minor species included black
cherry ( Prunus serotina)，eastern white pine ( Pinus
strobus)，white spruce (Picea glauca)，black spruce
(Picea mariana)，and other hardwoods． Tree diameter
at breast height ( DBH )，total height，crown length
( top to the base of crown)，and average crown width
were recorded for each tree ( Solomon et al．，2002 ) ．
The numbers of trees in each plot ranged from 12 to
109． Mean tree diameters were from 2. 2 to 19. 2 cm
and mean tree total heights were from 2. 7 to 17. 1
m (Tab． 1) ．
Tab． 1 Summary statistics of the stand
variables across the 50 plots
Stand variables Mean Std． Min． Max．
Number of trees per plot 34. 4 17. 7 12 109
Stand density /( tree·m － 2 ) 0. 61 0. 54 0. 18 2. 52
Stand basal area /(m2·hm － 2 ) 50. 2 14. 3 10. 0 79. 6
Stand mean DBH /cm 11. 5 3. 9 2. 2 19. 2
Stand mean height /m 11. 6 3. 7 2. 7 17. 1
Stand mean crown length /m 4. 4 1. 5 1. 6 8. 4
Stand mean crown width /m 1. 1 0. 3 0. 5 1. 9
Using the methods of refined nearest neighbor
statistic，Ripleys K-function (Ripely，1976) analysis
and pair correlation function analysis，24 plots were
classified as completely spatial random ( CSR ) point
pattern，17 plots as regular point pattern，and 9 plots as
clustered point pattern out of the 50 plots ( Li et al．，
2007) ． The classification was used as the basis for
modeling the spatial patterns in this study．
2. 2 Methods
2. 2. 1 Theoretical background of Gibbs model
(pairwise interaction process model) Suppose the
spatial point pattern in a mapped forest plot A has N
events ( trees) ． Denote X = {Xi = ( ui，vi)∈A ∶ i = 1，
2，…，N} as the set of coordinates ( ui，vi ) of the
trees． The forest plot A is a sampling window within a
much larger forest region． Therefore， the events
( trees) are the points of a partial realization of a
planar point process within A． For the patterns of N
trees in a bounded region A，the pairwise interaction
processes ( i． e．，pair-potential processes) are suitable
stochastic models． The joint density for these processes
generally takes the form of
f( x) = β
N
CN! {exp －1≤ i ＜ j≤N Φ( Xi － Xj }) ，
(1)
where · is the Euclidean distance between trees i
and j，β is a parameter determining the intensity of the
process，Φ(·) is a pair potential function，and C is a
normalizing constant to make f( x) a density depending
on β and Φ(·) ．
The homogeneous Gibbs point process is a special
case of pairwise interaction processes． It is
characterized by the potential function Φ (·) and
another parameter α = － log(β) ． Although the explicit
form of the distribution is hard to obtain，the following
relationship holds for a homogeneous Gibbs point
process (Diggle et al．，1994)
λE0［Z(X)］ =
{E Z(X) [exp － α － Σ∞
i = 1
Φ( Xi ] }) ， (2)
where λ is the density of events，E0 is the expectation
with respect to the Palm distribution of the process，E
is the expectation，and Z(X) is any random function
of the process where X includes all events of the
process in the whole plane． For a strict definition of
Palm distribution，see Stoyan et al． (1987)，Tomppo
(1986) and Cressie (1993) ．
However，when it is conditional on the number of
points N，Gibbs distribution can be expressed as
f( x) = Z －1 {exp －
1≤ i ＜ j≤N
Φ( Xi － Xj }) ，(3)
where Z is the normalizing constant in the form of
Z = ∫AN [exp －1≤ i ＜ j≤N Φ( Xi － Xj ]) dX1，…，dXN ．
The pair potential function Φ(·) can be interpreted as
followings: Φ ( d ) ＞ 0 indicates that the probability
density for inter-event distance d is lower than that
under Poisson process and，thus，there is repulsion at
distance d; Φ ( d) ＜ 0 indicates that the probability
density for inter-event distance d is higher than that
211
第 5 期 金星姬等: 树木位置空间模式建模与预测
under Poisson process． Thus， there is attraction at
distance d; and Φ(d) = 0 indicates a Poisson process．
2. 2. 2 Pair potential functions The pair potential
function Φ(·) can be parameterized in various forms．
In this study，three parametric models were adopted．
The first one is called Very-Soft-Core ( VSC ) pair
potential function (Ogata et al．，1984):
Φ( r) = － lg{1 － exp［－ ( r / δ) 2］}， (4)
where δ is a scaling parameter． For the VSC pair
potential function，the range of interaction is infinite
for all distance r． The second one is called Diggles
pair potential function that has the following form
(Diggle et al．，1994):
Φ( r) =
－ lg(1 －［1 － ( r / δ) 2］2)， if r ≤ δ;
0， if r ＞ δ{ ．
(5)
where the scaling parameter δ defines the range of
interaction． The third one is called Ogata and
Tanemuras pair potential function ( Ogata et al．，
1985):
Φ( r) = － lg(1 +［τ( r / δ) － 1］
exp［－ ( r / δ) 2］)，τ ≥ 0，δ ＞ 0， (6)
where the scaling parameter δ defines the range of
interaction． Notice that equation(4) is a special case of
equation(6) when the parameter τ = 0． For equation
(5)，Φ ( r) $0 when r = δ representing the Poisson
process． For equation (6 )，Φ ( r) $0 when τ = δ / r，
which represents the Poisson process．
2. 2. 3 Parameter estimation method One of the
methods for estimating the parameters of the pair
potential functions is the Takacs-Fiksel method (Diggle
et al．， 1994 ) ． It is based on the fundamental
relationship in equation ( 2 ) ． Suppose θ is the
parameter vector of a pair potential function ［θ = δ for
the potential functions(4) and(5) or θ = ( τ，δ) for the
potential function(6)］． If a series of test functions Zk
( k = 1，2，…，m) are chosen to estimate both sides of
equation(2) ( denoting the left side as L^ k ( α，θ) and
the right side as R^ k ( α，θ ))， the sum of squared
differences of both sides can be minimized to estimate
the model parameters (α，θ) as，
S(α，θ) = Σ
m
k = 1
［L^ k(α，θ) － R^k(α，θ)］
2， (7)
where m is some arbitrary number no less than the
length of parameter vector (α，θ) ．
There is no agreed rule for choosing the function
Zk ． However，if the test function Zk is chosen to be
equal to Nj ( rk ) (Diggle et al．，1994)，which is the
number of events i with Xi － Xj ≤ rk，where Xj( j =
1，2，…，L) is the jth randomly chosen testing point，
the left side of equation (2) can be estimated as
L^ k(α，θ) = λ
2K( rk)， (8)
where K ( rk ) is Ripleys K-function ( Ripley，1976;
1977 ) ． The right side of equation ( 2 ) can be
estimated as
R^ k(α，θ) =
1
LΣ
L
j = 1
Nj( rk)
[exp － α － ΣN
i = 1
Φ( Xi － Xj ]) ． (9)
The choices of L，m and rk are arbitrary although
different values can be tried and compared． For the
data used in this study， we set the constants as
follows: 1) L = 2N because we wanted to ensure that
there were enough random points surrounding a tree，
and 2 ) m = 30 and rk = 0. 1k ( which made the
maximum value 3. 0 m) because we wanted there to be
enough but not too many trees within the plot which are
considered interacting with each other．
The optimization of equation ( 7 ) was performed
using a direct search polytope algorithm ( one available
subroutine is the UMPOL in the IMSL Math Library of
Fortran PowerStation 4. 0 ℃ 1994 － 95 Microsoft
Corporation) ． Details of the polytope algorithm can be
found in the literature (Nelder et al．，1965; Gill et al．，
1981) ．
2. 2. 4 Edge corrections Typically， the study
regions would be the sampled plots within a forest
stand． Each of these plots covers only a part of the
entire spatial pattern． Interactions can only be
observed among the trees within a plot， while
interactions are unobserved between the trees within
the plot and the trees outside the plots． It is known as
the edge or boundary effect． If no correction is taken
for the boundary effect in the point pattern analysis，
the estimation of the interactions is biased． In the
previous study ( Li et al．，2007 )， the toroidal edge
correction method was proven to be simple and
satisfactory compared with other available edge
correction methods such as buffer zone， border
311
林 业 科 学 49 卷
method，and weighting method (Moeur，1993; Haase，
1995) since the plots in this study were all square or
rectangular in shape，and were small in size． Thus，
the toroidal edge correction method was used in
estimating the parameters of the pair potential functions
Φ(·) ．
2. 2. 5 Simulation of Gibbs models-Metropolis
algorithm Given a pair potential function，the Gibbs
point process models can be used to simulate an
equilibrium point pattern． One of the methods was
originally developed by Metropolis et al． (1953)，and
modified by Wood (1968) ． The algorithm is a Markov
Chain Monte Carlo (MCMC) procedure and has been
applied in many studies ( Ogata et al．，1981; 1984;
1985; Diggle et al．，1994)，which is briefly described
in Appendix．
2. 2. 6 Goodness-of-fit test To evaluate the
goodness-of-fit for each model，500 realizations of the
Gibbs models corresponding to each fitted pair potential
function were performed． Ripleys K-function
( converted to L as L^ ( d ) = K^(d) /槡 π － d ) was
calculated for each realization，where d is the given
distance of magnitude between paired events． The
95% confident envelopes were constructed and
plotted． Then，the observed L statistics from the data
were compared with the simulated 95% confident
envelopes． If model fitting to the data was good，the
observed L statistics should lie within the simulated
95% confident envelopes (Ripley，1976) ．
2. 2. 7 Predictions of the parameters of Gibbs
point process models The Gibbs point process
models are fitted based upon detailed mapped data and
individual tree measurements of diameter and height．
In forestry practice， however， it is difficult and
expensive to obtain these kinds of data． Therefore，it is
desirable to develop empirical regression models for
predicting the parameters of the Gibbs models using the
available stand variables as predictors． Since most of
the forest variables are highly correlated ( e． g．，tree
DBH and height)，the multicolinearity problem may
exist among the predictor variables． In this study，the
principal component ( PC ) regression was used to
overcome the multicolinearity problem ( Neter et al．，
1990) ．
3 Results
3. 1 Parameter estimation and goodness-of-fit
tests for three pair potential functions
The parameters of three pair potential functions，
namely VSC［equation(4)］，Diggles［equation(5)］
and Ogata and Tanemuras ［equation ( 6 )］， were
estimated by the Takacs-Fiksel method (Diggle et al．，
1994) for the 50 plots． The results of the goodness-of-
fit test indicated: overall，82% － 84% of the 50 plots
were satisfactorily fitted by the Gibbs models with three
pair potential functions ( Tab． 2 ) ． The Diggles pair
potential function performed best for the CSR plots．
The VSC pair potential function performed best for the
regular plots． The Ogata and Tanemuras pair potential
function was slightly better than the other two for fitting
the clustered plots (Tab． 2) ．
Tab． 2 Summary of goodness-of-fit tests of spatial point patterns for the 50 plots
Spatial point pattern Number of plots
Unaccepted
VSC Diggle Ogata ＆ Tanemura
Plot % Plot % Plot %
CSR 24 3 12. 5 0 0. 0 1 4. 2
Cluster 9 4 44. 4 4 44. 4 3 33. 3
Regular 17 2 11. 8 4 23. 5 5 29. 4
Total 50 9 18. 0 8 16. 0 9 18. 0
3. 2 Developing prediction models for the
parameters of Gibbs model
The Gibbs model with the Diggles pair potential
function ［equation ( 5)］ was selected as the spatial
point process model since it had slightly better
goodness-of-fit results among the three pair potential
functions used in this study．
The Pearson correlation coefficients between the
dependent variables ( i． e．，the estimated parameters δ^
and α^ of the Gibbs model with the Diggles pair
potential function) and the available stand variables
were computed and tested for H0: ρ i = 0． The
411
第 5 期 金星姬等: 树木位置空间模式建模与预测
correlation coefficients and the P-values of the tests
were shown in Tab． 3． It was evident that both δ^ and α^
had significant correlations ( P-value ＜ 0. 05 ) with
most of the stand variables．
The ordinary least squares regression was first
used to regress the two parameters ( δ^ and α^ ) on the
available stand variables such as stand density，basal
area，mean diameter，mean tree height，mean crown
length，and mean crown width． It was found that
multicolinearity existed among the predictor variables，
especially among the mean tree diameter，mean tree
height， and mean of crown length ( results not
shown) ． Thus， principal component regression was
applied to overcome the multicolinearity problem
(Neter et al．，1990) ． The regression coefficients based
on the principal component regressions for predicting
the two parameters ( δ^ and α^ ) of the Gibbs model on
the available stand variables were list in Tab． 4． The
model R2 was 0. 57 for the δ^ model and 0. 65 for the α^
model． Therefore，the Gibbs model parameters can be
predicted using the two empirical regression models
when the values of the six stand variables are available
( either observed or estimated) ．
Tab． 3 Correlations among the two estimated parameters of Gibbs point process model with
the Diggles pair potential function and the available stand variables
δ^ α^
Density /
( tree·m － 2 )
Basal area /
(m － 2·hm － 2 )
Mean
DBH /cm
Mean
height /m
Mean crown
length /m
Mean crown
width /m
δ^ 1 0. 338 8 － 0. 640 9 0. 334 1 0. 719 4 0. 705 9 0. 504 7 0. 273 0
(P-value) ( － ) (0. 016 1) (0. 000 0) (0. 017 7) (0. 000 0) (0. 000 0) (0. 000 2) (0. 060 5)
α^ 1 － 0. 788 2 0. 358 6 0. 725 2 0. 730 1 0. 576 7 0. 338 6
(P-value) ( － ) 0. 000 0 0. 010 6 0. 000 0 0. 000 0 0. 000 0 0. 018 6
Density 1 － 0. 541 5 － 0. 875 9 － 0. 891 8 － 0. 651 2 － 0. 473 2
(P-value) ( － ) 0. 000 0 0. 000 0 0. 000 0 0. 000 0 0. 000 7
Basal area 1 0. 698 7 0. 660 1 0. 598 5 0. 197 3
(P-value) ( － ) 0. 000 0 0. 000 0 0. 000 0 0. 178 9
Mean DBH 1 0. 968 6 0. 734 1 0. 407 1
(P-value) ( － ) 0. 000 0 0. 000 0 0. 004 1
Mean height 1 0. 681 4 0. 379 4
(P-value) ( － ) 0. 000 0 0. 007 8
Mean crown length 1 0. 236 5
(P-value) ( － ) 0. 105 5
Mean crown width 1
(P-value) ( － )
Tab． 4 Coefficients of empirical regression models for predicting the two parameters ( δ^ and α^) of
Gibbs point process model with the Diggles pair potential function using the available stand variables
Dependent
variable
Intercept
Density /
( tree·m － 2 )
Basal area /
(m2·hm － 2 )
Mean
DBH /cm
Mean
height /m
Mean crown
length /m
Mean crown
width /m
Fit statistics
R2 P-value
Diggles δ^ － 4. 215 7 0. 286 2 － 0. 038 2 0. 188 6 0. 220 6 0. 067 9 － 0. 117 4 0. 572 ＜ 0. 000 1
Diggles α^ 0. 920 8 － 0. 786 2 － 0. 009 4 0. 015 0 0. 027 9 0. 064 0 － 0. 086 2 0. 645 ＜ 0. 000 1
3. 3 Prediction of the spatial point patterns of the
50 plots
Based on the predicted parameters ( δ^ and α^) of
the Gibbs model with the Diggles pair potential
function，the spatial point patterns of a plot can be
simulated using the Metropolis algorithm． To evaluate
the performance of the two empirical regression models
for predicting the two parameters δ^ and α^，the Ripleys
K-function was used to compare the observed data and
model prediction． First，the observed K-function was
computed using the tree data for each of the 50 plots．
Then，for each plot，a 95% confident envelope was
constructed from 200 simulations using the predicted
parameters of the Gibbs model obtained from the two
empirical regression models in Tab． 4． If the observed
K-functions are within the 95% confidence envelope，
the prediction of the Gibbs model parameters can be
claimed satisfactory． The results indicated that 100%
511
林 业 科 学 49 卷
of the CSR plots，71% of the regular plots，and 56%
of the clustered plots were predicted well by the two
empirical regression models． In total，about 81% of
the plots were predicted satisfactorily by these models．
3. 4 Three example plots
One plot from each of the three spatial point
patterns was arbitrarily selected as the example to show
the application of the prediction models for the Gibbs
point process model． The estimated and predicted two
parameters ( δ^ and α^ ) as well as the stand variables of
the three plots were shown in Tab． 5． The estimated
two parameters were obtained from fitting the Gibbs
model with the Diggles pair potential function to the
tree data of each plot． The predicted two parameters
were obtained from the empirical regression models in
Tab． 4 using the observed values of the six stand
variables of each plot． Then，the predicted parameters
δ^ and α^ were used in the simulation process． Fig． 1，2，
and 3 illustrate 1 ) the observed versus simulated
Ripleys K-functions， and 2 ) the observed versus
simulated tree locations for the three example plots，
respectively． The 95% confidence envelopes for the
simulated K-function were based on 200 simulation
runs，while the simulated K-function curves shown in
the figures were randomly chosen from the 200
simulations．
Tab． 5 Estimated (Est． ) and predicted (Pred． ) parameters of Gibbs point process model
with the Diggles pair potential function for the three example plots
Plot
Spatial point pattern
05
Clustering
26
CSR
57
Regular
Stand variables
Density 0. 720 0 0. 402 4 0. 755 1
Basal area 46. 895 2 62. 072 8 57. 645 7
Mean DBH 8. 553 5 13. 035 0 9. 140 5
Mean HT 9. 607 4 12. 911 0 9. 526 8
Mean crown length 3. 128 8 5. 268 5 2. 930 8
Mean crown width 1. 114 4 1. 389 0 1. 194 4
Diggles pair
potential function
Parameter Est． Pred． Est． Pred． Est． Pred．
δ^ － 2. 783 5 － 1. 985 8 － 1. 349 8 － 0. 969 3 － 3. 652 2 － 2. 316 0
α^ 0. 212 3 0. 413 2 0. 990 9 0. 792 4 0. 338 4 0. 271 3
Plot 05 represents a clustered spatial point
pattern． The simulated Ripleys K-function only
partially follows the observed K-function for plot 05．
On the other hand，the simulated K-functions follow
the observed K-functions quite well for both plot 26
( representing a CSR spatial pattern ) and plot 57
( representing a regular spatial pattern)，respectively．
Nevertheless， the simulated tree locations look very
similar to the observed ones for the three example
plots． The overall patterns are reproduced successfully
for the CSR， regular and clustered spatial patterns
(Fig． 1，2，and 3) ．
Fig． 1 Observed versus simulated spatial point patterns of plot 05 ( representing a clustering spatial point pattern)
611
第 5 期 金星姬等: 树木位置空间模式建模与预测
Fig． 2 Observed versus simulated spatial point patterns of plot 26 ( representing a CSR spatial point pattern)
Fig． 3 Observed versus simulated spatial point patterns of plot 57 ( representing a regular spatial point pattern)
4 Discussion
In general，the Gibbs point process model with
three pair potential functions used in this study
performed satisfactorily in modeling the spatial point
patterns of trees within the 50 plots． The Diggles pair
potential function had slightly better goodness-of-fit
results than the other two． However，for each spatial
point pattern，the three models performed differently．
The CSR and regular plots were modeled better than
the clustered plots by all three pair potential functions．
This finding is consistent with the findings in other
studies using different spatial models ( Stoyan et al．，
2000) ． However，the shortcomings of this study were
1 ) the size of the available spruce-fir plots was
relatively small，2) the number of plots for each of the
three spatial patterns ( CSR，regular and clustering)
was limited，and 3) there were no independent data or
plots available for validating the empirical regression
models for predicting the two parameters of the Gibbs
model． Therefore，more and large forest plots for each
point pattern，especially the clustered plots，may be
needed to further confirm and validate the results in
this study．
For estimating the parameters of the three pair
potential functions， the Takacs-Fiksel method was
employed with arbitrary choices of the factors［such as
L，m and rk in equation (8)］． These factors may be
further evaluated by comparing the simulated spatial
point patterns for the estimated pair potential functions
from different combinations of these factor values． The
Metropolis algorithm was utilized for simulating the
spatial point patterns by the Gibbs point process
models． Again，several factors ［Δ，X ( 0 )，and T］
were arbitrarily determined and may affect the
convergence to equilibrium and even the final
equilibrium state． These factors were not formally
compared with different choices in this study．
However，the choice of T = 1 000 seems sufficient for
these plots given the values of other factors (Δ = 0. 5 m
711
林 业 科 学 49 卷
and the initial state as a random pattern of a Poisson
process) ．
In this study， the Gibbs point process models
were applied to all spatial distributions of trees． It has
the advantage of being flexible to fit various spatial
point patterns． This is extremely useful if the spatial
pattern of the data is not certain． Because of this
general property， its performance varies among
different spatial patterns． For example，in this study，
it did not fit the clustered pattern as well as it fitted the
CSR or regular patterns． In addition， the empirical
regression models for predicting the two parameters of
the Gibbs models worked best for the CSR pattern，
followed by the regular pattern， and worst for the
clustered pattern． This was consistent with the
goodness-of-fit of the Gibbs point process models for
the spatial patterns． Therefore，it would be interesting
to compare the model fitting and prediction，especially
for the clustered plots，using the Gibbs point process
model against other spatial models such as Cox process
model or Poisson cluster process model ( Matern，
1971; Cressie，1993) in a future study．
5 Conclusion
The Gibbs point process model with three pair
potential functions were used to model the spatial point
patterns of trees in the 50 spruce-fir plots in the
Northeast，USA． The results indicated that the three
pair potential functions fitted 82% － 84% of the plots
well． The Diggles pair potential function performed
slightly better than the other two pair potential
functions ( i． e．，VSC and Ogata and Tanemuras) ． In
general，the plots of CSR and regular patterns were
modeled better than the clustered plots by all three pair
potential functions． Further， empirical regression
models were developed to predict the two parameters
( δ^ and α^ ) of the Gibbs point process model with the
Diggles pair potential function using the available
stand variables as predictors ( i． e．， stand density，
basal area，mean diameter，mean tree height，mean
crown length，and mean crown width) ． The simulation
results showed that 81% of the 50 plots were predicted
satisfactorily by these models，in which 100% of the
CSR plots，71% of the regular plots，and 56% of the
clustered plots were predicted satisfactorily by the
Gibbs models． Three example plots were selected to
illustrate that the patterns of the simulated tree
locations were very similar to the observed ones．
The spatial structures of a forest stand include its
horizontal and vertical structures． The horizontal
structure can be characterized by the relative locations
of trees over space and the distributions of tree
diameters． The vertical structure can be represented by
the height distribution and diameter-height
relationships． The spatial models developed in this
study can provide useful information on the horizontal
structures of forest stands，whereas the information on
the vertical structures of forest stands can be obtained
from a bivariate distribution model of tree diameter and
height． The combination of these models will offer
much more information on the three dimensional profile
and dynamics of forest stands to forest managers and
researchers．
参 考 文 献
Cressie N A C． 1993． Statistics for spatial data: revised edition． John
Wiley ＆ Sons，New York，900．
Diggle P J． 1983． Statistical analysis of spatial point patterns． Academic
Press，New York，148．
Diggle P J，Fiksel T，Grabarnik P，et al． 1994． On parameter estimation
for pairwise interaction point process． International Statistical
Review，62(1) : 99 － 117．
Dovciak M，Frelich L E，Reich P B． 2001． Discordance in spatial
patterns of white pine (Pinus strobus) size-classes in a patchy near-
boreal forest． Journal of Ecology，89(2) : 280 － 291．
Fox J C，Ades P K，Bi H． 2001． Stochastic structure and individual-tree
growth models． Forest Ecology and Management，154(1 /2) : 261 －
276．
Friedman S K， Reich P B， Frelich L E． 2001． Multiple scale
composition and spatial distribution patterns of the north-eastern
Minnesota presettlement forest． Journal of Ecology， 89 ( 4 ) :
538 － 554．
Gill P E，Murray W，Wright M． 1981． Practical optimization． Academic
Press，New York．
Haase P． 1995． Spatial pattern analysis in ecology based on Ripleys K-
function: Introduction and methods of edge correction． Journal of
Vegetation Science，6(4) : 575 － 582．
Hara T． 1992． Effects of the mode of competition on stationary size
distribution in plant populations． Annals of Botany， 69 ( 6 ) :
509 － 513．
Jensen J L，Moller J． 1991． Pseudolikelihood for exponential family
models of spatial point processes． Annals of Applied Probability，1
(3) : 445 － 461．
Kenkel N C． 1988． Pattern of self-thinning in jack pine: testing the
random mortality hypothesis． Ecology，69(4) : 1017 － 1024．
811
第 5 期 金星姬等: 树木位置空间模式建模与预测
Kleinschmidt S，Baskerville G L． 1980． Foliage weight distribution in
the upper crown of balsam fir． USDA Forest Service Research Paper
NE-455，8．
Kuuluvainen T， Pukkala T． 1987． Effect of crown shape and tree
distribution on the spatial distribution of shade． Agricultural and
Forest Meteorology，40(3) : 215 － 231．
Landau L D，Lifshitz E M． 1980． Statistical Physics: Part I． 3rd ed．
Revised and enlarged by Lifshitz E M and Pitaevskii L P． Translated
from Russian by Sykes J B，Kearsley M J． Butterworth-Heinemann，
Boston．
Li F，Zhang L． 2007． Comparison of point pattern analysis methods for
classifying spatial distributions of spruce-fir stands in the Northeast，
USA． Forestry，80(3) : 337 － 349．
Matern B． 1971． Doubly stochastic Poisson processes in the plane∥Patil
G P，Pielou E C，Waters W E． Statistical ecology I． Pennsylvania
State University Press，College Town，PA，195 － 213．
Metropolis N，Rosenbluth A W，Teller M N，et al． 1953． Equation of
state calculations by fast computing machines． Journal of Chemical
Physics，21(6) : 1087 － 1092．
Miller T E，Weiner J． 1989． Local density variation may mimic effects of
asymmetric competition on plant size variability． Ecology，70 (4) :
1188 － 1191．
Moeur M． 1993． Characterizing spatial patterns of trees using stem-
mapped data． Forest Science，39(4) : 756 － 775．
Nelder J A， Mead R． 1965． A simplex algorithm for function
minimization． Computer Journal，7(4) : 308 － 313．
Neter J，Kutner M H，Nachtsheim C J，et al． 1990． Applied linear
regression models． 3 rd ed． IRWIN，720．
Ogata Y，Tanemura M． 1981． Estimation of interaction potentials of
spatial point patterns through the maximum likelihood procedure．
Annals of the Institute of Statistical Mathematics Series B，33(2) :
315 － 338．
Ogata Y， Tanemura M． 1984． Likelihood analysis of spatial point
patterns． Journal of Royal Statistical Society Series B，46 ( 3 ) :
496 － 518．
Ogata Y，Tanemura M． 1985． Estimation of interaction potentials of
marked spatial point patterns through the maximum likelihood
method． Biometrics，41(2) : 421 － 433．
Pacala S W，Deutschman D H． 1995． Details that matter: the spatial
distribution of individual trees maintains forest ecosystem function．
OIKOS，74(3) : 357 － 365．
Penttinen A，Stoyan D，Henttonen H M． 1992． Marked point processes
in forest statistics． Forest Science，38(4) : 806 － 824．
Pretzsch H． 1997． Analysis and modeling of spatial stand structure:
Methodological considerations based on mixed beech-larch stands in
Lower Saxony． Forest Ecology and Management， 97 ( 3 ) :
237 － 253．
Ripley B D． 1976． The second-order analysis of stationary point
processes． Journal of Applied Probability，13(2) : 255 － 266．
Ripley B D． 1977． Modeling spatial patterns (with discussions) ． Journal
of Royal Statistical Society Series B，39(2) : 172 － 192．
Rouvinen S，Kuuluvainen T． 1997． Structure and asymmetry of tree
crowns in relation to local competition in a natural mature Scots pine
forest． Canadian Journal of Forest Research，27(6) : 890 － 902．
Srkk A，Tomppo E． 1998． Modeling interaction between trees by
means of field observations． Forest Ecology and Management，108
(1 /2) : 57 － 62．
Solomon D S，Zhang L． 2002． Maximum size-density relationships for
mixed softwoods in the northeastern USA． Forest Ecology and
Management，155(1 /3) : 163 － 170．
Stoyan D，Kendall W S，Mecke J． 1987． Stochastic geometry and its
applications． John Wiley ＆ Sons，New York．
Stoyan D， Penttinen A． 2000． Recent applications of point process
methods in forestry statistics． Statistical Science，15(1) : 61 － 78．
Stoyan D，Stoyan H． 1998． Non-homogeneous Gibbs process models for
forestry - a case study． Biometrical Journal，40(5) : 521 － 531．
Tomppo E． 1986． Models and methods for analysing spatial patterns of
trees． Seloste: Malleja ja menetelmiae puiden tilajaerjestyksen
analysoimiseksi． Communicationes Instituti Forestalis Fenniae，138:
65．
Wood W W． 1968． Monte Carlo studies of simple liquid models∥
Temperly H N，Rowlinson J S，Rushbrooke G S． Physics of Simple
Liquids． Amsterdam，North-Holland，115 － 230．
Appendix: Simulation procedure of Gibbs models -Metropolis algorithm
Suppose there are N trees on a plot A，assuming to have periodic boundaries ( i． e．，square or rectangular) ．
These trees interact with each other according to a pair potential function，Φ(·) ． Assume the state of the N-tree
system is X( t) = {Xi( t)，i = 1，2，…，N} at stage t，where xi( t) is the ith trees location on A． The total potential
energy at the stage t is
Un(X( t)) = －
1≤ i ＜ j≤N
Φ( X ( t) i － X ( t) j ) ．
Then，a trial state X( t) = {Xi( t)，i = 1，2，…，N} is chosen such that a randomly selected tree X
*
r ( t) is
uniformly distributed within Xr( t) ± Δ． It defines a small square centered at Xr( t) with side length Δ on A，while
all other trees keep their locations on A the same as at the state t． The total potential energy for the trial state is
Un(X( t)) = －
1≤ i ＜ j≤N
Φ( X ( t) i － X ( t) j ) ．
Compare Un(X( t)) with Un(X( t)) as the following to decide the next state at stage ( t + 1): ( i) if Un(X( t))≤
Un(X( t))，the next state X( t + 1) = X( t)，and ( ii) if Un (X( t)) ＞ Un (X( t)) and if ξ≤exp［Un (X( t)) －
911
林 业 科 学 49 卷
Un(X( t))］，where ξ is a random number uniformly distributed on (0，1)，the next state is X( t + 1) = X( t);
Otherwise，the next state is X( t + 1) = X( t) ． For a sufficiently large T，the state X( T) will be a realization of
equilibrium state under the Gibbs distribution characterized by the pair potential function Φ(·) ． The choice of Δ，
initial state，and T affect the rate of convergence to the equilibrium． In this study，Δ was chosen to be 0. 5 m，the
initial state was a set of randomly generated locations of N trees on the plot A equivalent to a Poisson pattern，and
the maximum number of iterations T was set to 1 000．
(责任编辑 石红青)
021