阅读说明：本技术 一种建立大区域条件下单木生长模型的方法 (Method for establishing single-tree growth model under large-area condition ) 是由 吴兆飞 张春雨 王娟 范春雨 赵秀海 于 2019-10-30 设计创作，主要内容包括：本发明公开了一种建立大区域条件下单木生长模型的方法,包括如下步骤：数据收集、计算相关驱动因子、变量筛选、确定模型形式、分析和结论。本发明取得的有益效果为：(1)将气候因子包含在单木生长模型中,解决了气候差异对单木生长的影响；(2)单木生长模型具有较好的预测性,预测精度高。(The invention discloses a method for establishing a single-wood growth model under a large-area condition, which comprises the following steps: data collection, calculation of relevant drive factors, variable screening, determination of model form, analysis and conclusions. The invention has the following beneficial effects: (1) the climate factor is contained in the single-wood growth model, so that the influence of climate difference on the growth of the single wood is solved; (2) the single-wood growth model has good predictability and high prediction precision.)

一种建立大区域条件下单木生长模型的方法

技术领域

本发明涉及林业技术领域，更具体的说是涉及一种建立大区域条件下单木生长模型的方法。

背景技术

单木生长模型是基于影响树木生长的主要驱动因子定量描述单株树木生长过程的模型，在林业中具有广泛的应用。当前单木生长模型研究中最关键的一个问题就是如何确定树木生长驱动因子，从而更好地对树木生长状况和潜力进行预测。树木生长受到立木自身大小、林分条件、竞争因子和立地状况等因素的影响，并且作用强度与树种有关。

目前，由于研究条件的限制，大都集中在局部小范围的调查，缺少大区域调查数据(指的是不仅仅局限于某一个小范围，而是在整个地区尺度上建立模型)，大区域尺度条件下的单木生长模型仍然较少，气候因子未能包含在单木生长模型中，无法解决区域气候差异对单木生长的影响，这很大程度上降低了模型的应用范围，进而限制了重要树种的可持续经营和管理。

因此，如何提供一种建立单木生长模型的方法以解决区域环境对单木生长的影响是本领域技术人员亟需解决的问题。

发明内容

有鉴于此，本发明提供了一种建立大区域条件下单木生长模型的方法。

为了实现上述目的，本发明采用如下技术方案：

一种建立大区域条件下单木生长模型的方法，包括如下步骤：

(1)数据收集；

(2)计算相关驱动因子：

所述相关驱动因子包括树木自身大小、竞争指标和立地条件；

(3)变量筛选：

根据VIF＝1/(1-R^2)计算方差膨胀因子VIF，其中，R²指可决系数r²的算数平均值，即指用某一个相关驱动因子做因变量，对其它所有变量做回归以后得到的可决系数r²的算数平均值；

式中，y_i为因变量的真实值，为因变量的平均值，

为因变量的预测值；

VIF用来检验各驱动因子之间是否存在多重共线性：当0＜VIF＜10时不存在多重共线性，当10≤VIF＜100时存在较强的多重共线性；当VIF≥100时存在严重多重共线性，仅保留不存在多重共线性的相关驱动因子；

(4)确定模型形式：

确定模型为：ln(VI)＝b+f₁(SIZE)+f₂(COMP)+f₃(SITE)，

其中，VI为每年蓄积增量，SIZE为树木自身大小，COMP为竞争指标，SITE为立地条件，包括土壤条件和气候环境；

每年蓄积增量为单株树木每年的蓄积增长量(m³/year)。

VI代表生产力，在以往的研究中通常利用断面积增量表征生产力，本发明采用每年蓄积增量来表征，能更加全面的反映树木生长。

(5)分析：

利用逐步回归进行相关驱动因子的筛选和模型模拟，利用决定系数R²、平均误差

和均方根误差RMSE对模型进行评价，R²越大，

和RMSE越小则表明模型建立越精确；

逐步回归法是指将模型中逐个引入变量，每次引入对Y影响显著的变量，且每引入一个变量,对老变量逐个检验，将变为不显著的变量逐个剔除，并且将最终得到AIC信息统计量最小的模型作为最优模型。在R语言软件中利用step()函数实现。

(6)结论：

根据建模结果对单木生长模型进行分析，得出单木生长模型拟合关系式。优选的，所述步骤(1)包括如下步骤：

(11)确定研究区域，根据研究区域植被类型、物种组成、森林面积、气候条件、群落受干扰、破碎化程度以及经纬度坐标差异，均匀布设采样调查点，要求调查点数量与研究区域的实际情况相匹配；

(12)在所述调查点选取中心点，建立以中心点为圆心、半径为17.85m的样圆，并记录中心点的经纬度坐标、海拔、坡度和坡向；

(13)从正北方向开始，按顺时针方向测量所述样圆内乔木物种的胸径DBH，将DBH≥5cm的乔木记为目标树，记录所述目标树的物种名称、相对位置、DBH、高度和冠幅；

(14)对目标树进行年轮条取样、编号保存，待年轮条干燥后进行固定并打磨至可看清树木年轮为止，测定目标树近5年的生长量，结合胸径树高式和立木材积公式，求得每年蓄积增量。

胸径树高式是指在树高测量不便的情况下，利用胸径树高式进行计算，即：根据包含胸径的公式来计算树木高度。

立木材积式是指根据树高和胸径计算树木蓄积量的公式。

优选的，所述步骤(14)中年轮条取样的具体操作为：使用内径为5.15mm的生长锥对目标树进行年轮条取样，在树高1.3m处由北向南进钻，进钻方向与树干保持垂直，避开巴结或其它异常树皮处，进钻深度应超过髓心2～3cm，取出年轮条。

优选的，所述步骤(2)中树木自身大小包括胸径DBH和高度H；竞争指标包括林分密度指数SDI、胸高断面积密度BA、大于目标树的胸高断面积和BAL、Hegyi指数及混交度M_j；气候条件包括年均温度MAT、年均降水量MAP、最冷月温度CMT、生长季平均温度GST和总降水GSP；土壤状况包括总氮TOTN和碳氮比CNrt。

SDI是将现实林分的林木株数换算到标准平均直径时所具有的株数，是反映林分直径和株数密度的综合指标。BAL和Hegyi指数是与距离有关的两个竞争指标，BAL指的是目标树周围5m范围内相邻木中大于目标树个体的胸高断面积和，由于大植株相对于小植株有一定的竞争优势，大树胸高断面积可用于反映目标树所处的竞争环境中的非对称竞争。Hegyi指数考虑了目标树与周围相邻木之间的距离，利用相邻木的胸径和它们之间的距离来表示。混交度M_j指的是距离目标树最近的4株个体树种，异种个体所占的比例，用以反映目标树周围的混交状况。

优选的，所述步骤(2)中计算方法如下：

式中：SDI为林分密度指数，N为现实林分每公顷株数，D₀为标准木平均直径，D为现实林分的平均直径；

N为每公顷森林的树木株数；D₀为标准木平均直径，选取方式不一，我国一般为10cm，在本研究中由于起测径阶为5cm，所以选择15cm为标准木平均直径；D为现实林分的平均直径，为调查的所有目标树胸径的平均值。

胸高断面积指树木距地面以上1.3米处树干的横断面积。

式中：BA为林分胸高断面积密度，指每公顷林分所有树木胸高断面积和(cm²·ha^-1)，n指样圆内树木的株数，ba_i为第i株树木的胸高断面积，S指样圆的面积；

BAL为目标树周围5m范围内相邻木中大于目标树个体的胸高断面积之和；

相邻木是指目标树周围5m范围内的所有个体。

大于目标树个体指胸径大于目标树即可。

胸高断面积指树木距地面以上1.3米处树干的横断面积。

式中：HEGYI为目标树j的Hegyi指数，ba_i和ba_j为竞争木和目标树的胸高断面积，d_ij指的是目标树j和竞争木i的距离，n表示5m范围内竞争木的株数；

竞争木是指胸径大于目标树的相邻木。

式中：M_j指的是目标树j的混交度，当相邻木i与目标树j为相同的树种时，v_ij的值记为0，当相邻木i与目标树j为不同的树种时，v_ij的值记为1。

根据调查点的经纬度坐标跨度，利用KNMI Climate Explorer获取覆盖整个研究区的温度和降水数据，利用ArcGIS中的反距离插值法提取各调查点的年均温度MAT、年均降水量MAP、最冷月温度CMT、生长季平均温度GST和总降水GSP。

经由上述的技术方案可知，与现有技术相比，本发明取得的有益效果为：(1)将气候因子包含在单木生长模型中，解决了气候差异对单木生长的影响；(2)单木生长模型具有较好的预测性，预测精度高。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。

图1附图为本发明提供的建立单木生长模型的步骤流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

时间：2017年

(1)数据收集：

(11)确定研究地点，研究地点位于我国东北主要温带森林分布区，包括内蒙古自治区、黑龙江省、吉林省和辽宁省的8座山脉，由北至南依次为：大兴安岭、小兴安岭、完达山、老爷岭、张广才岭、长白山、哈达岭和龙岗山，地理位置为119°48'～134°01'E，39°42'～53°22'N。该地区气候类型复杂多样，主要为温带季风气候，冬季时间长，1月平均温度在-20℃以下，降雨主要集中于夏季，森林植被主要分布在山脉和丘陵。根据各山脉植被类型、物种组成、森林面积、气候条件、群落受干扰以及破碎化程度等差异，采用机械布点的方法，将8个温带森林分布区划分成了446个调查点。

(12)在每个调查点，选取中心点架设罗盘仪，调平并确定正北方向，建立以中心点为圆心、半径为17.85m的样圆。使用手持GPS和罗盘仪确定实际调查点的经纬度坐标、海拔、坡度、坡向等地形因子。

(13)使用胸径尺、卷尺、测高器和测距仪等工具，从正北方向开始，按顺时针方向测量所述样圆内乔木物种的胸径DBH，将DBH≥5cm的乔木记为目标树，记录目标树的物种名称、相对位置(北偏角，距离)、DBH、高度和冠幅。

(14)使用内径为5.15mm的生长锥对目标树进行年轮条取样，在树高1.3m处由北向南进钻，进钻方向与树干保持垂直，避开巴结或其它异常树皮处，进钻深度应超过髓心2-3cm；对年轮条进行编号保存，待所述年轮条干燥后进行固定并打磨至可看清树木年轮为止，并测定目标树近5年的生长量。根据目标树2017年的胸径以及近5年的生长量，计算出各目标树2012年的胸径，结合各地区各树种的胸径树高式和立木材积公式，求得每年蓄积增量。

公式可以根据《中国立木材积表》(2017)、《吉林省立木材积、出材率表》(2015)、《黑龙江省立木材积表》(1981)进行查阅。

(2)计算相关驱动因子：

相关驱动因子包括树木自身大小、竞争指标和立地条件。

树木自身大小包括胸径DBH和高度H；竞争指标包括林分密度指数SDI、胸高断面积密度BA、大于目标树的胸高断面积和BAL、Hegyi指数及混交度M_j；气候条件包括年均温度MAT、年均降水量MAP、最冷月温度CMT、生长季平均温度GST和总降水GSP；土壤状况包括总氮TOTN和碳氮比CNrt。

式中：SDI为林分密度指数，N为现实林分每公顷株数，D₀为标准木平均直径，D为现实林分的平均直径；

BAL为目标树周围5m范围内相邻木中大于目标树个体的胸高断面积之和；

式中：BA为林分胸高断面积密度，指每公顷林分所有树木胸高断面积和(cm²·ha^-1)，n指样圆内树木的株数，ba_i为第i株树木的胸高断面积，S指样圆的面积；

式中：HEGYI为目标树j的Hegyi指数，ba_i和ba_j为竞争木和目标树的胸高断面积，d_ij指的是目标树j和竞争木i的距离，n表示5m范围内竞争木的株数；

式中：M_j指的是目标树j的混交度，当相邻木i与目标树j为相同的树种时，v_ij的值记为0，当相邻木i与目标树j为不同的树种时，v_ij的值记为1。

根据调查点的经纬度坐标跨度，利用KNMI Climate Explorer获取覆盖整个研究区的温度和降水数据，利用ArcGIS中的反距离插值法提取各调查点的年均温度MAT、年均降水量MAP、最冷月温度CMT、生长季平均温度GST和总降水GSP。

(3)变量筛选：

根据VIF＝1/(1-R^2)计算方差膨胀因子VIF，其中，R²指可决系数r²的算数平均值，即指用某一个变量做因变量，对其它所有变量做回归以后得到的可决系数r²的算数平均值。

式中，y_i为因变量的真实值，为因变量的平均值，

为因变量的预测值。

VIF可以用来检验各驱动因子之间是否存在多重共线性：一般认为，当0＜VIF＜10时不存在多重共线性，当10≤VIF＜100时存在较强的多重共线性；当VIF≥100时存在严重多重共线性，对具有较强多重共线性的相关驱动因子进行剔除。

根据各变量的检验结果，发现林分密度指数SDI、生长季平均温度GST、总降水GSP、年均温度MAT、年均降水量MAP、最冷月温度CMT以及林分胸高断面积BA的VIF值大于10。通过依次删除VIF值较大的变量，最终保留下10个影响因子。具体结果见表1。

表1变量筛选结果

(4)确定模型形式：

确定模型的基本形式为ln(VI)＝b+f₁(SIZE)+f₂(COMP)+f₃(SITE)，式中，VI为每年蓄积增量，SIZE为树木自身大小，COMP为竞争指标，SITE为立地条件，包括土壤条件和气候环境。

(5)计算机编程实现分析过程：

利用R软件中的逐步回归进行相关驱动因子的筛选和模型模拟，利用决定系数R²、平均误差

和均方根误差RMSE对模型进行评价，R²越大，

和RMSE越小则表明模型越好。

(6)分析结果并得出结论：

根据逐步回归结果对蒙古栎单木生长模型进行分析，得到拟合关系式：ln(VI)＝-5.2198+0.0653*DBH+0.0630*H-0.0007*SDI-0.1853*BAL-0.0094*Hegyi+0.0396*Mj+0.0011*GSP+0.0187*GST-0.0527*TOTN+0.0151*CNrt

其中，决定系数R²为0.5951，平均误差E为-0.0047，均方根误差为0.0315，这表明所建立的蒙古栎单木生长模型具有较好的预测优度；RMSE值为0.0315，这表明包含区域环境差异因子的单木生长模型具有较高的精度。

在建立的包含区域环境差异因子的单木生长模型中，选取的相关驱动因子基本都包含在模型当中，这表明这些相关驱动因子对树木生长有显著影响。胸径和高度是树木生长的关键相关驱动因子，对蒙古栎生产力的影响都是显著正向的。林分密度指数、大于目标树的胸高断面积和、Hegyi指数三个竞争指标对蒙古栎单木生产力的影响均为显著负向，随着竞争的加剧，单木生长受到限制，生产力降低。此外，混交度对蒙古栎生长有促进作用，这说明混交度越大，不同树种由于生态位差异会提高资源利用率，促进树木生长。生长季温度和降水对蒙古栎都有显著影响，生长季温度和降水的升高会促进蒙古栎的生长。总氮对蒙古栎的影响为负向的，碳氮比对蒙古栎的影响为正向。

因此，在进行蒙古栎经营和管理的过程中，为了提高蒙古栎的生产力，应在允许的条件下尽可能降低蒙古栎周围的竞争，提高混交状况，并根据实际的气候条件进行预测经营。

对比例1

步骤(2)中不考虑气候条件进行生长模型分析，其余步骤同实施例1。

所得拟合关系式为：ln(VI)＝-4.3817+0.0621*DBH+0.0711*H-0.0007*SDI-0.1739*BAL-0.0105*Hegyi+0.0538*M_j

其中，决定系数R²为0.5798，平均误差E为-0.0048，RMSE值为0.0355，实施例1的RMSE值相比对比例降低了11.3％，这表明包含区域环境差异因子的单木生长模型具有较高的精度。

实施例2

采用本发明技术方案建立了兴安落叶松、黑桦、紫椴、色木槭、辽东栎的单木生长模型，具体评价结果见表2。由表2结果可知：无论对于何种树种，包含区域环境差异因子模型的R²值都有不同程度的提升，同时，RMSE值也有不同程度的降低，说明包含区域环境差异因子的模型有更好的预测优度。

表2各树种模型评价结果

各树种拟合关系式如下：

含区域环境差异因子：

兴安落叶松

ln(VI)＝-6.5830+0.0563*DBH+0.0004*SDI-0.2338*BAL+0.3024*M_j-0.00004*Hegyi+0.0352*CNrt+0.0002*GSP+0.0833*GST

黑桦

ln(VI)＝-7.370+0.0516*DBH+0.0423*H17+0.0002*SDI-0.0283*BAL-0.0002*Hegyi-0.0349*TOTN+0.0710*CNrt+0.0011*GSP+0.0652*GST

紫椴

ln(VI)＝-6.8430+0.0407*DBH+0.0736*H17-0.0195*BAL-0.00001*Hegyi-0.0235*TOTN+0.0435CNrt+0.0011GSP+0.0360*GST

色木槭

ln(VI)＝-6.314+0.0524*DBH+0.0900*H17-0.0011*SDI+0.1752*Mj-0.00001*Hegyi+0.0472*CNrt-0.0001*GSP+0.0234GST

辽东栎

ln(VI)＝-8.736+0.0876*DBH+0.1448*H17-0.1989*BAL-0.00001*Hegyi-0.0229*TOTN+0.0041*GSP+0.0417GST

不含区域环境差异因子：

兴安落叶松

ln(VI)＝-5.5017+0.0540*DBH+0.0005*SDI-0.2475*BAL+0.4913*M_j-0.00004*Hegyi+0.0588*CNrt

黑桦

ln(VI)＝-5.5762+0.0526*DBH+0.0463*H17+0.00002*SDI-0.0305*BAL-0.0001*Hegyi-0.0692*TOTN+0.0760*CNrt

紫椴

ln(VI)＝-5.8200+0.0406*DBH+0.0794*H17-0.1629*BAL-0.00001*Hegyi-0.0212*TOTN+0.0483*CNrt

色木槭

ln(VI)＝-5.8664+0.0532*DBH+0.0875*H17-0.0011*SDI+0.1595*Mj-0.00001*Hegyi+0.0497*CNrt

辽东栎

ln(VI)＝-5.4992+0.0782*DBH+0.1490*H17-0.3242*BAL-0.00001*Hegyi-02659*CNrt

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。

本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。