在植物生长过程中，随着植株的生长，越来越多空间被占据，导致可用的资源不足，最终引起一部分个体死亡，即自然稀疏(Westoby，1984； Hynynen，1993)。在林业上，自然稀疏法则常被用来构造相对密度指数(Drew et al.，1979)、密度控制图(Newton，2012)等，因此，分析自然稀疏规律对于林业管理者在实践上有很大的指导意义。

对于自然稀疏法则，讨论较多的是Reineke(1933)方程和Yoda等(1963)3/2法则。2位学者认为，对于一定的植物大小，必定存在一个最大的密度，即在一个林分生长演替过程中，存在着一条最大密度线。Reineke(1933)通过密度与平均直径的关系研究，发现最大密度线存在着固定斜率-1.605，而Yoda等(1963)研究发现，在自稀疏过程中，植物平均质量与密度之间的关系存在着固定斜率-1.5，即著名的3/2法则。之后，有研究者利用实地数据证明了Reineke和Yoda的结论： 最大密度线的斜率不随树种、林龄和立地质量的变化而变化，或者变化较小(Gadow，1986； Pretzsch et al.，2005)。然而，也有研究者对该结论提出了质疑，认为最大密度线的斜率是变化的(Bi et al.，1997；Pretzsch，2006)。Weiskittel等(2009)认为，最大密度线斜率估计方法不同是导致存在争议的主要原因之一。

估计最大密度线的方法主要有： 1)人为根据数据点波动范围的上限手绘最大密度线，并计算出斜率(Yoda et al.，1963； Drew et al.，1977)；2)最小二乘法(Wilson et al.，1988； Yang et al.，2002)，该方法简单，用得也比较多； 3)分位数回归法(Cade et al.，1999)；4)主成分分析和压轴回归分析法(reduced major axis regression)(Wilson et al.，1999；Bégin et al.，2001)；5)平分线回归法(bisector regression approach)(Newton，2006)； 6)随机边界法(Bi，2001；2004)等。Zhang等(2005)对几种估计最大密度线的方法做了比较，发现随机边界法要优于最小二乘法、分位数回归法和主成分分析法。

在分析最大密度线时，取样往往来自不同的立地、经营措施(初始密度)等。而如果采用以上几种方法对这些数据进行分析，一般都先假设最大密度线的截距和斜率在不同立地条件或经营措施下一致； 但是此假设存在着一定的争议。混合效应模型分析方法能够很好地解决这个问题，即将立地环境或经营措施当作随机效应。VanderSchaaf等(2007)利用混合效应模型估计火炬松(Pinus taeda)最大密度线，发现混合效应模型要优于其他方法； 但是混合效应模型仅给出参数的点估计值，并未能体现出参数估计值的不确定性(Li et al.，2012)。

分层贝叶斯模型法(hierarchical Bayesian method)是分析分层数据和评价参数不确定性的一个很好方法； 而且贝叶斯法综合利用了先验信息和样本信息，先验信息是在进行统计推断时不可缺少的一个因素，可以来自历史文献资料或者主观信念，这在林业研究工作中是很重要的(张雄清等，2014)。杉木(Cunninghamia lanceolata)是重要的用材树种，分布于我国南方大部分地区，在我国南方林区的林业可持续发展中占有重要地位。国内一些学者对杉木最大密度线也做了些研究。吴承祯等(2000)在前人研究植物种群自然稀疏法则的基础上，根据生长的密度理论和若干假设，提出描述杉木种群密度变化的理论模型； 陈水松等(2006)研究了20年生杉木间伐试验林自然稀疏的变化，并构建了密度控制管理图； Sun等(2011)利用几种不同估计方法分析了杉木自然稀疏线(材积-密度)，并选择出了最优估计方法等。但是对于分层贝叶斯法在杉木最大密度线的研究未见报道。本文以福建杉木密度试验林为研究对象，利用分层贝叶斯法分析杉木最大密度线，以期为杉木自然稀疏规律研究提供一种新思路。

1 试验地概况及数据整理

试验地位于福建省邵武市卫闽林场，武夷山北段中山山脉东南侧山区(117°43′E，27°05′N)。主要地貌为低山、高丘，海拔一般在250～700 m之间。该地区属于亚热带季风气候，气温温和湿润，年均温17.7℃，1月和7月的平均温度分别为6.8 ℃和28 ℃。年降雨量1 768 mm，年平均相对湿度为82%，地理条件比较适宜杉木生长，是杉木的中心产区。

该试验林使用1 年生苗木于1982 年造林，采用随机区组试验设计，5 个密度：A 密度(2 m × 3 m)、B 密度(2 m × 1.5 m)、C 密度(2 m × 1 m)、D 密度(1 m × 1.5 m)、E 密度(1 m × 1 m)，每个密度3 次重复，共15 个样地，每个样地面积600 m²。1983—1988年间逐年调查，1988 年后隔年调查，直到2010年。样地内每株树进行编号，对样地内林木进行每木检尺。孙洪刚等(2010)研究认为，当样地内枯损株数未达到总株数的10%时，杉木林分未开始进行自然稀疏，因此在拟合最大密度线时将这些数据剔除。本研究也采用相同方法对数据进行筛选，最终A密度的3个样地和B密度的1个样地都被排除掉，剩下11个样地，组成样本量为76的数据集。杉木林分中各变量统计值见表 1。

2 研究方法 2.1 最大密度线方程

Reineke(1933)根据林分平方平均直径(Q)和每公顷株数(N)的关系，提出最大密度线方程：

(1)

式中： a，b为待估参数； b为最大密度线的斜率。

在Reineke(1933)方程中，认为b是一个常数，为-1.605。之后，Yoda等(1963)在此方程理论基础上，提出了植物平均生物量和密度的关系，并得到了著名的3/2法则。White(1981)认为树木的生物量和材积有着非常密切的关系，因此，在林业上一般将平均生物量用平均材积(V)代替，得到以下方程：

(2)

在本研究中，利用Reineke(1933)方程分析杉木最大密度线，原因有以下3点： 1)Q是描述林分平均大小的常用变量之一，而且易获取； 2)Q与N和林分断面积关系很密切； 3)林木材积一般通过材积表查询得到。然而通过材积表获取的材积存在一定的误差，如果利用材积分析最大密度线[式(2)]，那么就将此误差引进最大密度线方程中(Curtis et al.，2000)，估计不准确。

2.2 贝叶斯理论

贝叶斯法是基于贝叶斯定理而发展起来的一种系统的统计推断方法。贝叶斯推断的基本原理是将待估参数的先验分布与样本信息结合，再根据贝叶斯定理得出后验分布信息，之后根据此后验分布信息去推断未知参数的分布(张雄清等，2014)。令y=(y₁，y₂，y₃，…)为数据集，θ=(θ₁，θ₂，θ₃，…)为待估的参数向量，则根据贝叶斯理论，其基本公式为：

(3)

式中： p为概率分布函数或者密度函数。

在贝叶斯方法中，通过概率分布 p(θ|y) 描述参数θ的不确定性，进而估计参数θ。根据贝叶斯条件概率，θ的条件概率分布为：

(4)

贝叶斯统计中，在给定样本数据y下参数θ的条件分布 p(θ|y) 即所需求参数的后验分布， p(y|θ) 是在给定θ下y的似然函数，p(θ)是θ的先验分布。

在本研究中，根据最大密度线方程(1)，得到下式：

(5)

式中： g(lnN:a，b)=a+blnN 。所以方程(3)可以写成：

(6)

其中 p(样本|a，b) 为:

(7)

2.3 分层贝叶斯

初植密度对最大密度线是否有影响长期以来一直饱受争论。在本研究中，为分析初植密度对杉木最大密度线的影响，引入了分层贝叶斯模型。模型形式有3种： 1)在截距处考虑初植密度的随机效应，而斜率作为固定效应不变(M1)； 2)在斜率处考虑初植密度的随机效应，而截距作为固定效应不变(M2)； 3)同时在截距和斜率处考虑随机效应(M3)。3种方程形式如下：

(8)

(9)

(10)

式中： u_0i，u_1i为随机效应参数，而且u_0i～N(0，σ₀²)，u_1i～N(0，σ₁²)，ε～N(0，σ²)。对于分层贝叶斯模型，模型的总方差有2个来源：组内波动性(within-group variability)和组间波动性(between-group variability)，在本研究中，即不同初植密度间的波动性和相同初植密度内的波动性。

2.4 先验分布

贝叶斯方法的一个重要特点是参数在抽样前就具有先验分布(Gelman et al.，2004)。在上述杉木最大密度线方程中，首先要为参数a，b构造合适的先验分布。有些学者提出利用无信息先验分布，对于无信息先验分布，一般选择均值为0、方差足够大的正态分布，使其能够覆盖整个数据范围的波动情况(Ellison，2004)。当然也可以选择有信息先验分布作为贝叶斯方法中的先验分布，这些先验信息可以来自文献资料或者主观信念。本研究中，对于斜率b，通过查阅相关文献，搜集树木自然稀疏的最大密度线的斜率值(表 2)，构造先验信息分布。而对于截距a，则利用无信息先验分布。参数的先验信息分布见表 3。

3 模型评价

对于模型的评价，一般采用决定系数(R²)和均方根误差(RMSE)作为模型的拟合统计量指标。均方根误差公式为：

(11)

式中： y_i 为实际值；$\bar{y}$_i，y_i 分别为预测值和平均值；n为观察个数。

在本研究中，同时还利用DIC统计量判断贝叶斯模型拟合的效果(Spiegelhalter et al.，2002)。与RMSE一样，DIC越小，说明模型拟合越好。

(12)

式中： DBAR=E_θ{-2lg[p(y|θ)]}， 表示模型拟合数据的优劣； PD为模型中参数的有效个数， PD=DBAR+2lg[p(y|θ-)]， 表示模型的复杂度。

对于贝叶斯参数估计，利用WinBUGS软件(Spiegelhalter et al.，2003)。在该软件中，通过Gibbs抽样(Chib et al.，1995)，利用马尔科夫链蒙特卡洛方法完成估计。在本研究中通过R2WinBUGS(Sturtz et al.，2005)连接R软件和WinBUGS软件完成，同时输出图形。在估计参数时，为了保证迭代收敛和得到稳定的参数后验概率值，迭代次数设为10万次，并去掉前面的5 000次退火(burn-in)迭代。

4 结果与分析

在杉木最大密度线方程中，分别利用普通贝叶斯模型和分层贝叶斯模型估计最大密度线的截距和斜率，得到估计值和模型评价指标，见表 4。虽然分层贝叶斯模型的R²比普通贝叶斯模型大，均方根误差RMSE比普通贝叶斯模型的小，但是参数估计值a和b与普通贝叶斯相近；而且，也发现普通贝叶斯的DIC绝对值小于M1模型、大于M2模型。综合DIC，R²和RMSE比较可发现，M2模型要比其他模型好。

基于模型M2，得到图 1。杉木最大密度线的截距和斜率服从一定的正态分布，这也很好地解释了最大密度线的截距和斜率不是固定不变的，而是存在着一定的不确定性。此外，根据模型M2，得到杉木最大密度线图，见图 2。通过分层贝叶斯估计所拟合的最大密度线相对较好，但是对于最大密度线的上界有一定的低估。

另外，通过表 4也发现，根据T-检验，随机效应方差( σ₀²，σ₁² )不管是在截距上( σ₀² =0.008，SD=0.029)还是斜率上( σ₁² =0.003，SD=0.016)，其t值都小于1.96，差异不显著，由此可推断初植密度对杉木最大密度线的斜率和截距都没影响。因此，在表 4中未列出M3的值。在模型M1中可发现，总方差为0.035，其中初植密度间的波动性所带来的方差为0.008； 换而言之，只有22.86%的方差是由于不同初植密度所造成的，而剩下77.14%的方差是由于相同初植密度内重复测量所造成的。这同样的问题在模型M2中也可发现，由初植密度所带来的方差仅占总方差的10%。

5 结论与讨论

对于最大密度线争论较多的是斜率是否为常数。有些学者认为斜率是常数，不随树种、林龄、立地等因子的变化而变化(Puettmann et al.，1993； Tang et al.，1995； Yang et al.，2002)； 然而，有些学者却认为最大密度线的斜率随着初植密度的不同而异(Turnblom et al.，2000； VanderSchaaf et al.，2007)。本研究发现，初植密度并没有对杉木最大密度线的斜率造成影响，且初植密度在截距上的随机效应不显著，因此也可理解为初植密度对杉木最大密度线的截距影响不显著，与其他一些研究(Weiskittel et al.，2009； Leon，2011)不同。由于初植密度的随机效应不论在截距上还是在斜率上都不显著，因此最大密度线模型的不确定性在很大程度上是由相同初植密度内不同样地的重复测量所造成的。Zeide(2001)发现林龄和环境的变化能够使得树木的自然稀疏发生变化。Harrington(1986)和Weiskittel等(2009)研究发现赤桦木(Alnus rubra)的最大密度线对坡面、土壤水分、干燥度等很敏感。因此，最大密度线在不同的样地存在一定的不确定性，这也可以从本文的图 1反映出来。

尽管对于最大密度线的估计已经有很多种方法，但是这些方法都存在一定的局限性。如通过手绘的方法画最大密度线比较主观、武断(Drew et al.，1977)； 分位数回归(Scharf et al.，1998； Cade et al.，1999)估计最大密度线波动性比较大，甚至数据微小的变化都会使得估计结果产生较大的波动(Scharf et al.，1998； Cade et al.，2003)； 最小二乘法描述的是数据的平均状态，这显然与实际的情况不符(孙洪刚等，2010)。Bi等(2000)利用随机边界法分析最大密度线，得到较理想的结果。然而这些方法忽略了最大密度线的波动性，而且也没有考虑重复测量间的相关性。在本研究中，提出利用分层贝叶斯模型分析杉木最大密度线，不仅能够很好地反映出杉木最大密度线的不确定性，而且也考虑了重复测量产生的相关性。另外，利用分层贝叶斯模型，引入最大密度线的斜率的先验信息分布，估计结果更加可靠，这也是贝叶斯方法估计的一大特点。本文所估计的杉木最大密度线的斜率的后验概率平均值为-1.882(M2)，该值与WBE模型(West et al.，1997； 1999)计算所得的-2相近。而且，Mohler等(1978)研究10个树种纯林的最大密度线的斜率，得到这10个树种的斜率范围为-2.29～-1.750，本研究所得到的值属于这个范围。因此，利用分层贝叶斯模型分析杉木最大密度线能够为杉木最大密度线的研究提供一种新思路。当然，也可以在杉木人工林最大密度线方程中加入立地质量、环境变量等因素，使得先验分布更精确，进而提高杉木人工林最大密度线预测的可靠性。