林分结构中最基础、最重要的是林分直径结构规律。在林分内，各种大小直径树木按径阶的分配状态，称作林分直径结构(stand diameter structure)，亦称林分直径分布(stand diameter distribution)。林分直径结构无论在理论还是在实际上都是最基本的林分结构，不仅是因为林分直径便于测定，而且因为林分内各种大小直径树木的分配状态直接影响树木的树高、干形、材积、材种、树冠和生物量等因子的变化(吴可等，2010)。大量研究表明，上述各因子的结构规律与林分直径结构规律紧密相关，因此研究林分直径结构分布特征在理论上可为森林经营及测树制表方法和技术提供依据。

许多概率密度函数或模型被用来描述林分直径的分布特征，如Normal模型、Lognormal模型、Gamma模型、Johnson’s S_B模型、β模型和Weibull模型(Bailey et al., 1973；Burkhart et al., 1974；Hafley et al., 1977)。其中，Weibull三参数分布模型由于具有能灵活地拟合一系列曲线和偏度等级、相对容易估计出参数、封闭的累计密度公式等优点(Bailey et al., 1973；Bowling et al., 1989；Matney et al., 1982；Rennolls et al., 1985；Zarnoch et al., 1991), 常用来拟合和预测林分直径结构分布变化。

被广泛用来估计Weibull三参数分布模型的方法有最大似然法、矩法、百分位法和最小二乘回归法(Lei, 2008)，而预测Weibull三参数分布模型的方法主要有：1) 参数预测法(PPM)，直接拟合概率密度函数的参数；2) 参数回收法(PRM)，通过直径分布的矩来求解概率密度函数的参数；3) 参数百分位法(PCT)，用树木直径百分位的方法求解概率密度函数的参数。预测Weibull分布三参数的参数回收法就是估计Weibull分布的矩法，而参数百分位法就是以林木直径的百分位数与Weibull分布函数参数的数学关系式求解分布参数，它们都可以用来估计和预测Weibull函数参数及参数变化，以展示林分直径结构的动态变化趋势。

长期以来，这些方法分别被用来估计和预测不同树种和森林类型的林分直径分布，如采用参数预测方法研究的森林类型和树种有短叶松(Pinus banksiana)人工林(Smalley et al., 1974)、火炬松(Pinus taeda)(Clutter et al., 1984；Feduccia et al., 1979)、铁杉(Tsuga chinensis)和花旗松(Pseudotsuga menziesii)混交林(Little，1983)、挪威云杉(Picea abies)人工林(Kilkki et al., 1989；Siipilehto，1999)等，采用参数回收法研究的森林类型和树种有火炬松(Cao et al., 1982；Matney et al., 1982)、长叶松(Pinus palustris)(Leduc et al., 2001)、天然火炬松林(Burk et al., 1984)、阿巴拉契亚山脉硬阔林(Matney et al., 1987)等，采用参数百分位法研究的森林类型和树种有长叶松人工林(Lohrey et al., 1976)、火炬松人工林(Baldwin et al., 1987)、火炬松(Knowe, 1992)等。Liu等(2004)比较了参数预测法、参数回收法和参数百分位法预测黑云杉(Picea mariana)人工林直径分布的优劣，但没有比较不同参数估计方法拟合直径分布的效果。我国采用参数预测法预测林分直径分布的研究较多，如北京栎类(张雄清等，2009)、蒙古栎(Quercus mongolica)(陈新美等，2008)、油松(Pinus tabulaeformis)人工林(孟宪宇，1985)、杜仲(Eucommia ulmoides)人工林(李荣伟等，2000)、天然阔叶林(陈昌雄等，2004)、日本落叶松(Larix kaempferi)(方精云等，1987)、杉木(Cunninghamia lanceolata)人工林(周国模等，1992)等；采用参数回收法和参数百分位法求解Weibull三参数的研究较少，仅见魏柏松等(1995)采用参数回收法求解湿地松(Pinus elliotii)人工林的Weibull三参数、孟宪宇等(1991)采用参数回收法预测长白落叶松(Larix olgensis)直径分布收获模型；而采用多种方法的比较研究则更少，仅方子兴(1993)利用杉木标准地资料比较了最大似然法、矩法、百分位法和回归法4种Weibull参数估计方法拟合直径分布的优劣，但是没有比较参数预测法、参数回收法和参数百分位法预测直径分布变化的情况。

蒙古栎又称柞树或柞木，在我国主要分布于东北和华北地区，是我国温带地区落叶阔叶林及针阔混交林的主要树种(于顺利等，2000)。蒙古栎用途广泛，具有很高的经济价值，而且还具有保持水土、涵养水源的作用。为了提高森林质量，科学开展蒙古栎次生林经营非常重要，而研究其林分直径结构分布特征，比较Weibull直径分布参数估计和预测方法在蒙古栎次生林经营中的适用性，了解不同方法在蒙古栎林分直径分布参数估计和参数预测的精确性，对蒙古栎林经营具有非常重要的参考价值。本文利用157块蒙古栎纯林样地资料，比较最大似然法、矩法和百分位法估计Weibull三参数分布函数的效果，并用参数预测法、参数回收法和参数百分位法预测Weibull三参数直径分布结果，利用Kolmogorov-Smirnov(K-S)检验、误差指数和一些常用的模型估计检验方法比较不同方法估计和预测蒙古栎林分直径分布的优劣，同时比较不同方法估计和预测蒙古栎林分直径分布参数的精确性，为更好开展蒙古栎林经营提供理论依据和技术参数。

1 数据来源与研究方法 1.1 数据来源

数据来源于2009年测定的吉林省蒙古栎林分固定样地数据(N=157块，树木株数为13 068)，每块样地面积为0.06 hm²。主要调查因子包括平均胸径、方位角、水平距、平均年龄、平均树高、优势高、郁闭度、地位级、坡位、坡向、坡度、海拔、土层厚度和土壤类型等。平均年龄、平均树高是每块样地抽取3株平均木进行测量后的平均值，平均胸径为胸高断面积平均直径, 优势高是每块样地中最高的6株树的平均值。表 1所示为蒙古栎样地林分和单木水平因子统计。

1.2 分布模型和参数估计以及分布变化预测模型方法

Weibull三参数分布的概率密度函数f(x)为：

$\begin{array}{c} f\left(x \right) = \left({\frac{\gamma }{\beta }} \right){\left({\frac{{x - \alpha }}{\beta }} \right)^{\gamma - 1}}\exp \left[ { - {{\left({\frac{{x - \alpha }}{\beta }} \right)}^\gamma }} \right], \\ \alpha \le x \le \infty, \beta > 0, \gamma > 0. \end{array}$

(1)

式中：x为林木实测胸径；α，β，γ分别为Weibull分布模型的位置、尺度和形状参数。

Weibull三参数估计的3种方法如下。

1.2.1 最大似然法(MLE)

最大似然法是Weibull分布模型参数估计的一种常用方法，估计精度较高，但是其对计算的要求也较高。在最大似然法估计过程中，如果位置参数α为负值，则本研究设置α值为样地最小观测胸径的0.5倍。Weibull分布的其他两参数可以设计为：z_i=x_i-a，其中x_i为第i株林木胸径，a为林分最小胸径。最大似然法公式为：

$\begin{array}{c} \frac{1}{\gamma } + \frac{1}{n}\sum\limits_{i = 1}^n {\ln {z_i} - \left[ {\left({\sum {z_i^\gamma \cdot \ln {x_i}} } \right)/} \right.} \\ \left. {\left({\sum {z_i^\gamma } } \right)} \right] = 0; \end{array}$

(2)

$\beta = {\left({\frac{1}{n}\sum\limits_{i = 1}^n {z_i^\gamma } } \right)^{\frac{1}{\gamma }}}.$

(3)

式中：n为林分中各样地树木的株数。

根据式(2)(方子兴，1993)迭代求出参数γ，代入式(3)(方子兴，1993)求出参数β。为了预测林分直径分布变化，将最大似然估计得到的样地Weibull分布3个参数分别建立与其林分因子的线性回归关系，因此可以通过实测或预测林分因子的变化直接预测分布参数的变化，这就是参数预测方法(PPM)。

1.2.2 矩法(MOM)

Weibull概率分布密度函数的参数与林木胸径的矩之间存在数学关系，矩法(如算数平均直径为一阶矩，均方直径为二阶矩)可以理解为林木直径的一阶矩和均方直径的二阶矩与Weibull函数参数的数学关系。对于Weibull三参数分布，采用林木胸径的前二阶矩来估计参数α，β，γ。矩法的主要步骤为：1) 建立回归方程预测林分最小直径、算术平均直径和均方直径；2) 用式α=kx_min决定α的取值，式中k的取值可以是0，1/3，1/2，2/3，1(Bailey et al., 1973；Knoebel et al., 1986)，本研究取1/2；3) 建立被预测分布的前二阶矩与算术平均直径和均方直径的公式(Liu et al., 2004)如下：

$\beta = \frac{{\bar x - \alpha }}{{\Gamma \left({1 + \frac{1}{\gamma }} \right)}};$

(4)

$\bar x_{\rm{q}}^2 - \alpha \left({2\bar x - \alpha } \right) - \frac{{\left({\bar x - \alpha } \right)}}{{\Gamma {{\left({1 + \frac{1}{\gamma }} \right)}^2}}}\Gamma \left({1 + \frac{2}{\gamma }} \right) = 0.$

(5)

式中：Γ是gamma函数；x_q是均方直径；x是平均直径。当k已知时，可求出α，则可以用二分法求出式(5)中的γ，α和γ已知后，再代入式(4)即可求出β。

1.2.3 百分位法(PM)

百分位法与矩法类似，即以林木胸径的不同百分位作为Weibull分布参数的函数(Borders et al., 1987；Knowe，1992)。本研究利用第0，25，50和95的直径分布的百分位来估计Weibull分布的三参数，方程(Liu et al., 2004)如下：

$\alpha = \frac{{{n^{0.3333}}{x_0} - {x_{50}}}}{{{n^{0.3333}} - 1}};$

(6)

$\gamma = \frac{{2.343088}}{{\ln \left({{x_{95}} - \alpha } \right) - \ln \left({{x_{25}} - \alpha } \right)}};$

(7)

$\begin{array}{c} \beta = \frac{{\alpha \Gamma \left({1 + \frac{1}{\gamma }} \right)}}{{\Gamma \left({1 + \frac{2}{\gamma }} \right)}} + {\rm{sqrt}}\left\{ {\frac{{{\alpha ^2}}}{{{\Gamma ^2}\left({1 + \frac{1}{\gamma }} \right)}}\left[ {{\Gamma ^2}\left({1 + \frac{1}{\gamma }} \right) - } \right.} \right.\\ \left. {\left. {\Gamma \left({1 + \frac{2}{\gamma }} \right)} \right] + \frac{{\bar x_{\rm{q}}^2}}{{\Gamma \left({1 + \frac{2}{\gamma }} \right)}}} \right\}. \end{array}$

(8)

式中：n为林分中各样地林木株数；x₀，x₂₅，x₅₀，x₉₅分别为直径分布的第0，25，50和95的百分位值；x_q为均方直径；Γ为gamma函数；当式(6)中α < 0时，设置α为0；sqrt为开方。

为了预测林分直径分布的变化，参数预测法是由最大似然估计获取的Weibull分布参数与林分因子(例如林分年龄和密度等)建立线性回归模型，而参数回收法和参数百分位法是通过建立矩法和百分位法估计的Weibull分布参数与林分因子的线性回归关系，然后依据样地不同时间林分调查因子(例如林分年龄和密度等)来预测Weibull分布三参数的变化或林分直径分布的动态变化，这就是广泛使用的参数预测法(PPM)、参数回收法(PRM)和参数百分位法(PCT)预测林分直径动态变化的方法。

1.3 模型检验和评价

本研究采用2个指标来评价估计Weibull三参数3种方法的优劣。

第1个指标是Kolmogorov-Smirnov(K-S)检验值，公式如下：

$D = \mathop {Sup}\limits_{ - \triangleleft x < \infty } \left| {{F^*}\left(x \right) - F\left(x \right)} \right|.$

(9)

式中：D为K-S统计值；Sup为最大值；F^*(x)林分中各样地直径分布的理论分布；O(x)林分中各样地直径分布的累积概率函数。

第2个指标是预测误差指数(e)，公式如下：

$e = \sum\limits_j^m {\left| {{P_j}\left(x \right) - {O_j}\left(x \right)} \right|.} $

(10)

式中：m为每块样地的径阶数；P_j(x)为林分中各样地第j径阶的预测株数；O_j(x)为林分中各样地第j径阶的实际株数。K-S统计值和预测误差指数越小，说明预测效果越好。

本研究中，3种预测方法的评价使用确定系数(R²)和估计参数共线性的方差膨胀因子(VIF)。所有参数估计和评价方法利用Matlab完成。为了验证3种方法估计Weibull分布模型参数的预测效果，本研究将157块样地随机分为112块样地用于拟合模型参数、45块样地用于模型检验。

2 结果与分析

对157块样地数据，采用矩法、百分位法和最大似然法分别拟合每块样地的树木直径分布，并通过Kolmogorov-Smirnov(K-S)检验(P=0.05)，结果如表 2所示。由表 2可知，最大似然法的接受率为96.18%，矩法和百分位法的接受率分别为89.81%和82.80%。3种估计方法中最大似然法接受率最高，百分位法最低，矩法居中。

表 3列出了157块样地矩法、百分位法、最大似然法估计的Weibull三参数的平均值、K-S统计值的平均值和误差指数的平均值。由表 3可知，3种方法估计位置参数的变动在0.049 0~7.892 6之间，尺度参数的变动在0.064 5~66.286 3之间，形状参数的变动在0.783 4~15.251 1之间。

从平均误差指数来看，最大似然法的平均误差指数最低，为24.614 3，其次是矩法，为29.027 8，最后是百分位法，为29.590 2。从K-S统计值来看，最大似然法的K-S统计平均值最低，为0.097 8，百分位法最高，为0.149 2，矩法居中，为0.110 5。无论是平均误差指数还是K-S统计值，均为最大似然法最低。从表 3还可以看出，最大似然法的平均误差指数变动最小，标准差只有9.860 0，百分位法变动最大，标准差为22.381 5，矩法居中，标准差为11.690 8。因此，最大似然法在估计Weibull三参数时具有一定优势。

为预测林分分布结构的动态变化，对112块样地建立的参数回收法(PRM)、参数预测法(PPM)和参数百分法(PCT)估计的分布参数与林分因子的平均年龄、平均树高、优势高和对数林分密度[ln(林分密度)]进行多元线性回归，估计结果如表 4所示。由表 4可知，决定系数(R²)最低的是参数预测法的位置参数α方程，仅0.101 2；参数回收法的均方直径方程决定系数最高，达到0.863 9；参数百分位法的x₉₅百分位方程决定系数居中，为0.851 5。3种参数预测方法中最好的是参数回收方法，其方差膨胀因子(VIF)在1.214 6~2.306 3之间，均小于10，说明不存在多重共线性。

为了评价3种模型K-S值和误差指数的统计学意义，本文进行了配对t检验。在进行配对t检验之前，采用K-S值和误差指数检验拟合数据和检验数据是否满足正态性，结果表明，无论是拟合数据和检验数据的误差指数还是其K-S值，都支持正态假设。

从表 5可以看出，通过配对t检验，对于误差指数，矩法和百分位法的平均数对比在P < 0.05时差异不显著，矩法和最大似然法的平均数对比在P < 0.05时差异显著，百分位法和最大似然法的平均数对比在P < 0.05时差异显著。由前述(表 3)可知最大似然法的误差指数平均数最低，这说明最大似然法在0.05显著性水平时更加精确地拟合了数据。由表 6可知，任意2种方法的K-S统计值的对比在α=0.05时均差异显著。根据表 3，最大似然法的K-S统计值的平均数最小，因此通过配对t检验可知，最大似然法的误差指数平均数在显著水平为0.05时显著小于其他2种方法。

图 1是从157块样地中随机选出3块样地的直径分布直方图与最大似然法、矩法和百分位法估计的模拟图。从图 1a可以看出，最大似然法、矩法和百分位法均较好地拟合了单峰分布，3种方法基本没有差别；从图 1b可以看出，3种方法均较好地拟合了右偏的分布；从图 1c可以看出，3种方法均没有较好地拟合多峰分布，错过了高峰和低谷。相对来说，百分位法在拟合多峰分布时稍好一些。

表 7所示为45块样地3种预测Weibull模型参数的结果。通过K-S检验可知，45块检验样地中，PRM，PCT和PPM分别有19，16和21块样地没有通过检验，通过率分别为57.78%，64.46%和53.33%，PCT的通过率最高，其次PRM，最后是PPM。

PRM，PCT和PPM的平均误差指数分别为32.664 6，30.202 9和31.521 0。从表 8可以看出，通过配对t检验，PCT与PRM和PPM的误差指数平均数对比在P < 0.1时差异显著；也就是说，PCT的平均误差指数在P < 0.1时小于PRM和PPM，PCT在P < 0.1时比PRM和PPM更加精确地预测了数据。PRM，PCT和PPM的K-S统计值分别为0.194 9，0.188 7，0.207 9。从表 9可以看出，PRM和PPM、PCT和PPM的K-S统计值平均数对比在P < 0.1时差异显著；也就是说，PCT的K-S统计值平均数在P < 0.1时显著小于PRM和PPM。综合可知，PCT在预测数据时其精度较其他2种方法高。

为了解不同方法在预测不同径阶时的情况，本研究先将误差指数分径阶进行计算，然后求45个样地的平均值。图 2所示为3种方法分径阶的误差指数平均值情况。由图 2可知，误差指数平均值在整体上有随径阶增大而减小的趋势。3种方法在小径阶时，误差指数平均值均较高，而在大径阶时，误差指数均比小径阶时降低。本研究与Liu等(2004)研究结果相同，主要原因可能是小径阶的林木株数较少。最大的误差指数平均值在最小径阶处，即第6径阶。从变动来说，PRM的变动最大，而PCT的变动最小。

3 结论

本研究旨在比较Weibull三参数密度函数的矩法、百分位法和最大似然估计方法，以及参数预测法(PPM)、参数回收法的矩法(PRM)和参数百分位法(PCT)估计和预测蒙古栎纯林直径分布的差别。结果表明，3种估计方法均较好地估计了蒙古栎纯林的直径分布，K-S检验的接受率为82.80%~96.18%，其中最大似然法的接受率最高。

采用配对t检验分析了3种方法估计Weibull分布的误差指数和K-S统计值的平均数的差异，结果表明最大似然法可更加精确地估计蒙古栎纯林的直径分布。

采用PPM建立了Weibull函数的三参数与林分因子--平均年龄、平均树高、优势高和林分密度之间的回归模型，同时建立了PRM和PCT的参数回收法用以得到PRM和PCT三参数的解释性参数，如最小直径、平均直径、均方直径和直径频率百分位(第0，25，50，95)之间的回归模型，相关系数从0.101 2~0.863 9不等。研究还将回归方程计算得出的各参数代入Weibull分布，以预测直径分布。通过K-S检验可知，PRM，PCT和PPM的接受率分别为57.78%，64.45%和53.33%。通过配对t检验可知，PCT在显著水平为0.1时比PRM和PPM更加精确地预测了蒙古栎的直径分布。

另外，虽然3种估计方法在拟合多峰分布时均效果不好，但在拟合单峰分布和偏向分布时均取得了较好的效果。