白河林业局红松分布数量预估模型1)

王 烁 李凤日 甄 贞

(东北林业大学，哈尔滨，150040)

红松是整个北半球中、高纬度地区分布面积最大的3 大五针松(红松、西伯利亚红松和美国白松)之一，并以树干通直高大、材质优良、寿命长、生产力高、种子营养丰富等独特的用途与价值而闻名于世界，是目前世界上珍贵而稀有的多用途树种之一。同时，红松林作为东北地区重要而珍贵的森林资源和东北大平原农业发展的重要屏障，对人类生存区域气候环境的维护和调节起着非常重要的生态作用。从1949年起，红松就担负着沉重的商品木材生产任务，超负荷的过量采伐以及不当的保护措施，导致红松资源日益减少。我国东北东部山地的大面积原始红松混交林，经过近一个世纪的破坏和集中开发，到目前为止，可采森林资源接近枯竭［1］。为了保护以及合理开发利用红松资源，研究其分布有重要意义。另外，在现阶段的红松林管护承包经营中，红松资源的准确调查和估计已成为急需解决的问题。逻辑斯蒂模型(Logistic)适用于对“是非值”因变量的分析，泊松模型(Poisson)和负二项分布(Negative binomial，NB)模型适用于对计数因变量的分析，三者均属于广义线性模型，在社会、经济、林业等许多领域都有广泛的应用［2－14］。本研究根据最小二乘法(Ordinary lease square，OLS)，分别建立Logistic、Poisson 和NB 模型来预估天然红松分布概率和分布株数，并对所建模型进行拟合效果评价及独立性检验。与传统的全林调查方法相比，利用模型对红松株数分布进行预测，可以节省大量人力，物力和财力，有很好的应用价值。

1 研究区概况

白河林业局位于吉林省延边朝鲜族自治州安图县的南部，地处长白山脚下。其地理坐标为东经127°53'～128°34'，北纬42°01'～42°48'。该局经营区南北最长89 km，东西最宽56 km，全局总面积190 470 hm2，整个区域呈由东南向西北走向的狭长条形［15－17］。

白河林业局地带性顶极群落为以针叶树为主的针阔混交林。天然林主要林分类型为慢生阔叶混交林、落叶松林、针阔混交林、白桦林、中生阔叶混交林。由于人类的不断砍伐以及不合理利用，在天然林面积和蓄积中，针叶林(包括红松林、云杉林、长白赤松林、落叶松林、臭松林和针叶混交林)的面积和蓄积只占22.3%和22.8%［18］。

2 研究方法

2.1 数据来源与处理

白河林业局1990年结合森林资源二类调查，在全局范围内按2 km×1 km 网格，共布设了1 030 块固定样地(样地面积0.06 hm2)，目的是监测全局的森林资源动态。本研究采用2000—2010年的固定样地复测数据，从1 030 块固定样地中，剔除非林地、无林地、疏林地等样地258 块样地，选择了落入有林地的772 块固定样地作为基础数据。固定样地在林业局的分布如图1所示。每块样地调查因子包括海拔高、坡度、植被盖度、坡向、林龄、平均胸径、每公顷断面积等共21 个调查因子。其中，对坡向进行了数量化处理，即:平坡和北坡替换为0°，东北坡替换为45°，东坡替换为90°，东南坡替换为135°，南坡替换为180°，西南坡替换为225°，西坡替换为270°，西北坡替换为315°。在772 块固定样地中，随机抽取了70%即540 块固定样地数据用于建立模型，而剩余的232 块固定样地数据用于模型检验。

图1 白河地区样地分布图

2.2 模型的研究方法

2.2.1 模型变量的选择

计算机常用的模型选择方法有:正向逐步选择法、反向逐步选择法和混合逐步选择法。其中混合逐步选择法最普遍［19］。混合逐步选择法(combined stepwise)是根据所设的显著性标准，分别将变量加入到模型中去或删除掉。本研究使用SAS 9.2 软件，通过采用混合逐步选择法，将统计显著水平控制在α=0.05，得到了一些影响红松分布的变量，包括海拔高、坡度、植被盖度、植被高度、林龄、优势木胸径、平均胸径、平均树高、郁闭度、株数、每公顷断面积、活立木蓄积量、生物量地上部分、全树生物量、碳地上部分、全树碳含量，但是其中有些变量并没有意义，而有些对红松分布有意义的变量并没有包含其中。因此，以红松生长条件为基础，通过混合逐步选择法和林业经验对其进行增加与删除，最终保留海拔高、坡度、植被盖度、坡向、林龄、平均胸径和每公顷断面积7 个回归模型变量，各变量的统计量见表1。

表1 模型变量的统计量

2.2.2 Logistic 模型

Logistic 回归模型用来预测红松是否存在，模型的形式如下:

式中:P 为红松的分布概率，β0～β7为回归系数，X1表示海拔高，X2表示坡度，X3表示植被盖度，X4表示坡向，X5表示林龄，X6表示平均胸径，X7表示每公顷断面积。

2.2.3 Poisson 模型

Poisson 回归模型用来预测红松的分布数量，模型的形式如下:

式中:μ 为红松分布株数，β0～β7为回归系数，X1表示海拔高，X2表示坡度，X3表示植被盖度，X4表示坡向，X5表示林龄，X6表示平均胸径，X7表示每公顷断面积。

2.2.4 负二项分布模型

NB 模型是泊松模型的一种扩展形式，它一般是用来解决因变量的不均匀分布(即过度散布)的问题。当条件方差超过条件平均值时，负二项分布模型通过散布参数调整计数变量分布的异质性。负二项分布模型与Poisson 模型具有相似的概率密度方程，如公式(3)所示，

式中:Γ 表示伽马函数，μi服从伽马分布，即μi～Γ(α，β)，α 为形状参数，β 为比例参数。

2.3 模型检验

2.3.1 模型拟合效果评价

Logistic 模型拟合效果评价:在Logistic 模型中，有以下检验指标可以用来评价模型的拟合效果。

(1)Cox and Snell R2值:

(2)似然比检验:似然比检验是零假设检验之一，零假设为β=0，用于检验模型的所有系数是否为零。为了在两个以两组不同参数中似然函数不同比例的值为基础的假设中做出一个决定。似然比检验指标定义如下，

(3)Hosmer－Lemeshow 拟合优度指标:是一种类似于皮尔森χ2 统计量的指标，先将观测数据按照预测概率做升序排列，再将数据分为大致相同的n个组，其公式如下:

式中:G 代表分组数，ng代表第g 组的样本个数为第g 组的平均预测概率，yg代表第g 组样地的观测值。该统计量服从自由度为n－2 的χ2分布，若χ2检验不显著，即LH值很小，说明预测值与观测值间没有显著差别，表示模型可以很好地拟合数据，相反，LH值很大，模型拟合效果不理想。

(4)ROC 曲线:是根据一系列不同的二分类方式(分解值或决定阈)，以真阳性率(灵敏度)为纵坐标，假阳性率(特异度)为横坐标绘制的曲线。最靠近左上角的ROC 曲线的点是错误最少的最好阈值，其假阳性和假阴性的总数最小，所以ROC 曲线越靠近左上角，曲线下方面积越大，实验的准确性就越高，模型拟合效果越好。

(5)混淆矩阵:是将观测事件分为事件发生或不发生的频数表，如表2所示，对角线上的值(即a和d)为模型正确预测的样本个数，该表可以用于检验Logistic 模型的预测准确性。

表2 Logistic 模型拟合结果的混淆矩阵

Poisson 模型拟合效果评价:对于Poisson 模型，本研究计算以下检验指标来评价模型的拟合效果。

(1)AIC 值:AIC 是衡量统计模型拟合优良性的一种标准，它建立在熵的概念基础上，可以权衡所估计模型的复杂度和此模型拟合数据的优良性［20］。AIC 在比较模型的时候使用，值越小，说明模型的拟合效果越好。

(2)偏差分析:通过分析观测值间的偏差来分析模型拟合优度的方法，其公式如下:

2.3.2 独立性检验

模型的独立性检验是采用未参与建模的独立样本检验数据来对所建模型的预测性能进行综合评价。通常用平均偏差、平均绝对偏差和模型预测精度等偏差统计量作为评价模型预测能力的指标。

播种两周后，在田间按小区取水稻叶片混样，并提取DNA。水稻基因组DNA的提取采用李进波等［5］改进的CTAB法。

平均偏差:

平均绝对偏差:

模型预测精度:

式中:yi为实际值，为模型预估值，n 为样本数，t0.05为student t 表中α为0.05的概率值，=模型检验时，ME、MAE值越小，P 值越大则模型的预估精度越高，模型预估效果越好。

3 结果与分析

根据选取的7 个变量，建立了Logistic 模型用于对红松分布概率进行预测，模型形式为:

式中:Z1=0.896 7+0.000 5X1－0.065 7X2+0.030 3X3+0.002 4X4－0.028 8X5+0.091 8X6－0.034 6X7。P 为红松的分布概率，X1表示海拔高，X2表示坡度，X3表示植被盖度，X4表示坡向，X5表示林龄，X6表示平均胸径，X7表示每公顷断面积。

同时，表3列出了对模型参数估计值的显著性检验结果。除海拔高外，其余变量的p 值均小于0.05(α)，说明这些变量都通过了显著性检验，大多数变量对红松分布概率具有显著影响。其中，植被盖度和林龄的p 值小于0.000 1，说明这两个变量对天然红松分布概率具有最为显著的影响。在Logistic 回归模型的系数估计值中，坡度，林龄，每公顷断面积的估计值为正，说明这些变量与天然红松分布概率是正相关的。其中，坡度的值最大(0.065 7)，说明在其它变量保持不变的情况下，坡度每增加1°，红松分布概率提高约0.07。海拔高，植被盖度，坡向，平均胸径的系数估计值为负，说明这些变量与红松的分布概率成负相关。例如，在其它变量保持不变的情况下，如果平均胸径增加1 m3，红松分布概率会减少0.09。

表3 Logistic 模型参数估计值及其显著性检验

根据选取的7 个变量，建立了Poisson 模型用于对红松分布株数进行预测，模型形式为:

表4列出了对模型参数估计值的显著性检验结果。坡度，植被盖度，林龄，平均胸径，每公顷断面积的p 值均小于0.05(α)，说明它们对预测天然红松的分布株数影响显著。在Poisson 模型的系数估计中，每公顷断面积的估计值最大且为正，说明在其它变量保持不变的情况下，每公顷断面积每增加2 m2，大约多出现1 棵红松。平均胸径的估计值为负，其绝对值最大，说明平均胸径越小，天然红松的分布株数越多。

由于红松株数的分布过度分散(overdispersion)，本研究建立了负二项分布模型来处理因变量不均匀分布的问题，模型形式为:

式中:Z3=－0.465 7－0.000 5X1+0.041 1X2－0.014 3X3－0.000 6X4+0.023 2X5－0.107 1X6+0.051 7X7。

表5列出了对模型参数估计值的显著性检验结果。与Poisson 模型相似，坡度、植被盖度、林龄、平均胸径、每公顷断面积的p 值均小于0.05(α)，说明它们对预测天然红松的分布株数影响显著。估计值最大且为正的仍是每公顷断面积(0.052)，估计值为负，其绝对值最大的为平均胸径(－0.107)。

表4 Poisson 模型参数估计值及其显著性检验

表5 负二项分布模型参数估计值及其显著性检验

模型拟合效果评价:表6列出了Logistic 模型的似然比和Hosmer－Lemeshow 拟合优度检验。在似然比检验中，检验结果为p 值小于0.000 1，因此拒绝零假设，即模型系数不为0。在Hosmer－Lemeshow 拟合优度检验中，p 值为0.62，接受零假设，因此模型拟合效果很好。图2为模型的ROC 曲线，曲线下方的面积为0.76(＞0.5)，模型拟合效果较好。

表6 Logistic 模型的拟合效果检验

根据540 块固定样地红松分布概率的预测值，计算出期望值是0.2，因此，将概率界限设为0.2，如果一个观测的预测事件概率值等于或大于0.2，表示预测事件发生，否则为预测事件不发生。表7为Logistic 模型的拟合结果的混淆矩阵，总正确率为67%，主要错误是将红松出现事件预测为未出现事件(约占63.3%)。Logistic 模型的Cox and Snell R2值为0.165 2，在Logistic 模型中属于可接受的范围。

图2 ROC 曲线

表7 Logistic 模型拟合结果混淆矩阵

对于红松分布株数的预测(泊松模型和负二项分布模型)，由于观测样本中红松株数为0 的样地较多，本文考虑了泊松模型和负二项分布模型并比较两种模型，避免因变量的不均匀分布给模型拟合带来的影响。通过比较两个模型的CAI值，发现Poisson 模型的CAI值为1 609.69，而NB 模型的CAI值为1 060.45，因此，NB 模型的拟合效果更好。在红松分布株数预测模型拟合度检验中(表8)，Poisson 模型的P 值为0，拒绝零假设(模型很好地拟合观测数据);NB 模型的P 值为1，接受零假设。因此，与Poisson 模型相比，NB 模型的偏差明显变小，且与观测数据的拟合度较好。在似然比检验中，零假设为Poisson 模型和NB 模型拟合效果相同，P 值为0，说明拒绝零假设，即Poisson 模型和NB 模型拟合效果有显著不同，与NB 模型拟合效果更好的结论相符。

独立性检验:在Logistic、Poisson 和NB 模型中，利用未参与建模的232 块固定样地数据对所建模型进行独立性检验和预测精度检验。表9对3 个模型的独立性检验结果表明，Poisson 模型的平均偏差值较小，Logistic 模型的平均绝对偏差值较小，Logistic 模型的模型预测精度最高(70.23%)。Poisson和NB 模型预测略低于Logistic 模型，且精度相似(69.29%和68.16%)，说明对红松分布概率的预测要比对红松分布株数的预测准确率略高。表10 显示了Logistic 模型预测结果的混淆矩阵，对红松分布概率的预测正确率为71.1%，略高于模型拟合时的正确率(67%)，错误率最高的仍然出现在将红松出现事件预测为未出现事件(约占62.2%)。

表8 红松分布株数预测模型拟合度检验

表9 Logistic，Poisson 和NB 模型独立性检验结果

表10 Logistic 模型预测结果混淆矩阵

4 结论

根据最小二乘法分别建立全局的逻辑斯蒂、泊松和负二项分布模型，可以有效地预测和分析大区域内红松的分布概率及分布株数。与泊松模型相比，负二项分布模型更适合于解决不均匀分布的问题，负二项分布模型拟合数据效果明显优于泊松模型，独立性检验效果与泊松模型相似，因此，可以用负二项分布模型预测红松分布株数。红松一般生长于海拔在150～1 800 m 的地带，由于范围比较广，海拔对于红松生长的影响并不大。在Logistic 模型中海拔高的估计值为负，而在Poisson 模型中海拔高估计值为正，说明随着海拔高的增加，红松分布概率减小，但一旦确定有红松分布，则海拔越高，红松分布株数就会越多。林龄对于红松分布的影响是很显著的，随着林分年龄的增加，红松分布概率也增加，而且分布株数也相应增多。平均胸径对于红松分布的影响也是很显著的，随着林木生长，其平均胸径增加，林木间的竞争也会逐渐加强，林木的数量将会逐渐减少，所以红松分布概率减少，而且红松分布株数也会相应降低。

5 讨论

本研究根据7 个林分因子及地形因子建立了逻辑斯蒂、泊松和负二项分布模型，初步解决了红松在白河地区的分布概率及分布株数的预估问题。在独立性检验中，Logistic 模型的预测精度(70.23%)略高于Poisson 模型(69.29%)和NB 模型的预测精度(68.16%)，虽然精度尚可接受，但还有提高的空间。模型预测精度较低可能是由于建模过程中将每个样地视为独立的个体，而忽略了变量的空间自相关性。在今后的研究中，需要将样本的空间自相关性考虑在建模的过程中，建立预测红松分布概率和分布株数的地理加权模型(Geographically Weighted Regression，GWR)。GWR 是Fortheringham 等人在总结前人局部回归分析和变参数研究的基础上，基于局部光滑的思想提出的［21］，该模型将数据的空间位置嵌入到回归参数中，利用局部加权最小二乘方法进行逐点参数估计，其中权是回归点所在的地理空间位置到其它各观测点的地理空间位置之间的距离的函数。通过各地理空间位置上参数估计值随地理空间位置的变化情况，可以非常直观的探测空间关系的非平稳性，目前已被应用于城市地理学、气象学、森林学等诸多科学领域，并且应用在工业，农业，经济等方面，均有显著效果［22－26］。因此，在建模过程中考虑变量的空间自相关性将成为当代林业的一个重要发展趋势。

［1］ 傅俊卿.东北天然红松林资源现状与保护经营对策［J］.东北林业大学学报，2009，37(2):75－78.

［2］ 檀根甲，胡官保.Logistic 模型在植病流行学中应用的有关问题的探讨［J］.安徽农业大学学报，1998，25(2):138－141.

［3］ 胡喜生，范海兰，宋萍，等.改进Logistic 模型在城市人口预测中的应用［J］.北京大学学报，2008，9(4):370－373.

［4］ 孙洪刚，张建国，段爱国，等.5 种Logistic 模型在模拟杉木人工林胸高断面积分布中的应用［J］.林业科学研究，2007，20(5):622－629.

［5］ 马彪.青海省人口城市化Logistic 模型及其应用［J］.甘肃科技，2008，24(5):41－43.

［6］ 阎慧臻.Logistic 模型在人口预测中的应用［J］.大连工业大学学报，2008，27(4):333－335.

［7］ 周元满，谢正生，刘素青，等.Logistic 模型在桉树生长过程估计中的应用［J］.南京林业大学学报:自然科学版，2004，28(6):107－110.

［8］ 郭福涛，胡海清，马志海，等.不同模型对拟合大兴安岭林火发生与气象因素关系的适用性［J］.应用生态学报，2010，21(1):159－164.

［9］ Zhen Zhen，Li Fengri，Liu Zhaogang，et al.Geographically local modeling of occurrence，count，and volume of downwood in Northeast China［J］.Applied Geography，2013，37(5):114－126.

［10］ 李嘉竹，刘贤赵，李宝江，等.基于Logistic 模型估算水资源污染经济损失研究［J］.自然资源学报，2009，24(9):1667－1675.

［11］ 王远飞，张超.Logistic 模型参数估计与我国城市化水平预测［J］.经济地理，1997，17(4):8－13.

［12］ 雷渊才，张雄清.长白落叶松林分进界模型的研究［J］.林业科学研究，2013，26(5):554－561.

［13］ 张雄清，雷渊才，雷相东，等.基于计数模型方法的林分枯损研究［J］.林业科学，2012，48(8):54－61.

［14］ 孙龙，尚喆超，胡海清，等.Poisson 回归模型和负二项回归模型在林火预测领域的应用［J］.林业科学，2012，48(5):126－129.

［15］ 郭红，李凤日，龚文峰，等.基于GIS 的白河林业局景观格局演变研究［J］.森林工程，2009，25(4):1－5.

［16］ 张英瑞，李珊珊，张媛，等.吉林省白河林业局高保护价值森林判定研究［J］.当代生态农业，2013(Z1):124－130.

［17］ Xu Dong，Dai Limin，Shao Guofan，et al.Forest fire risk zone mapping from satellite images and GIS for Baihe Forestry Bureau Jilin China［J］.Journal of Forestry Research，2005，16(3):169－174.

［18］ 胡增军，贾炜玮，赵宇彤，等.白河林业局天然混交林的枯损率［J］.东北林业大学学报，2011，39(10):25－27.

［19］ 王济川，郭志刚.Logistic 回归模型:方法与应用［M］.北京:高等教育出版社，2001.

［20］ 宋喜芳，李建平，胡希远，等.模型选择信息量准则AIC 及其在方差分析中的应用［J］.西北农林科技大学学报:自然科学版，2009，37(2):88－92.

［21］ 覃文忠.地理加权回归基本理论与应用研究［D］.上海:同济大学，2007.

［22］ 吴玉鸣，李建霞.基于地理加权回归模型的省域工业全要素生产率分析［J］.经济地理，2006，26(5):748－752.

［23］ 刘琼峰，李明德，段建南，等.农田土壤铅、镉含量影响因素地理加权回归模型分析［J］.农业工程学报，2013，29(3):225－234.

［24］ 续秋霞，薛红，冯文娟，等.基于混合地理加权回归模型的我国人均GDP 分析［J］.纺织高校基础科学学报，2011，24(3):445－454.

［25］ 苏芳林.基于地理加权回归模型的县域经济发展的空间因素分析［J］.学术论坛，2005，5(5):81－84.

［26］ 俞路.基于GWR 模型的长江三角洲区域经济增长主导因素研究［J］.工业技术经济，2011(8):27－32.