基于哑变量的木荷人工林树皮厚度模型研究

张伟志 吴宏炜 黄光灿 田意 庄崇洋 江希钿

摘  要：在木荷人工林中，基于啞变量模型法构建树皮厚度模型，以提高预测模型参数的稳定性，为木荷人工林出材率计算和经济价值评估方面提供参考依据。基于福建省南平市213块木荷人工林固定样地的调查数据，选取19个包含胸高处、任意高度处、相对树皮厚度和去皮直径的理论方程作为构建木荷人工林树皮厚度的基础模型，使用R软件进行模型拟合分析，运用决定系数（R2）、均方根误差（Erms）、和方差（Ess）、赤池信息准则（AIC）、贝叶斯准则（BIC）等模型评价指标，结合含熵值的TOPSIS法筛选出拟合度较高的4个基础模型（M1、M10、M12、M16），进一步构建含龄组和立地质量哑变量的木荷人工林树皮厚度模型。用于建模的胸高处样本数为130个，任意高度处、相对树皮厚度和去皮直径样本数为2386个，用于检验模型拟合效果相对应的样本数分别为55个和1013个。结果表明：含龄组哑变量胸高处、任意高度处、相对树皮厚度模型（M20、M22、M25）的R2分别为0.9769、0.9214、0.9111，比基础模型（M1、M10、M12）的R2（分别为0.6981、0.5540、0.5056）提高了39.9%、66.3%、80.2%，有了明显提升;含龄组哑变量去皮直径模型M27的AIC为21.62，BIC为76.44，比其基础模型M16分别降低了95.8%、85.46%。通过含熵值的TOPSIS法对模型指标综合评价，并经模型配对t检验，这4个含龄组哑变量模型（M20、M22、M25、M27）拟合效果更佳，适合于福建木荷人工林树皮厚度预测。

关键词：树皮厚度;哑变量;龄组;立地质量;木荷

中图分类号：S758;S79      文献标识码：A

Abstract： In the artificial forest of Schima superba， a bark thickness model was constructed based on the dummy variable model method， in order to improve the stability of the prediction model parameters， and provide reference for the calculation of the timber yield and the evaluation of the economic value of the artificial forest of S. superba. Based on the survey data of 213 sample plots of S. superba plantation in Nanping City， Fujian Province， 19 theoretical equations including breast height， any height， relative bark thickness and peeling diameter were selected as the basic models for building the bark thickness of S. superba plantation. Based on R software model fitting analysis， determination coefficient （R2）， root mean square error （Erms）， the sum of squares due to error （Ess）， akaike information criterion （AIC） and bayesian information criterion （BIC） models as the evaluation index， and the entropy TOPSIS method to screen high fitting degree of four basic models （M1， M10， M12， M16） model， an age group and the site quality of S. superba plantation bark thickness of dummy variable model was further built. The number of samples at the breast height used for modeling was 130， the number of samples at any height， the relative bark thickness， and the peeling diameter was 2386， and the corresponding numbers for testing the model fitting effect was 55 and 1013 respectively. The results showed that R2 of the model with age dummy variable at breast height， any height and relative bark thickness （M20， M22 and M25） was 0.9769， 0.9214 and 0.9111， respectively， which was 39.9%， 66.3% and 80.2% higher than that of the basic model （M1， M10 and M12）， and AIC and BIC of M27 with age dummy variable was 21.62 and 76.44， respectively， it was 95.8% and 85.46% lower than that of the basic model M16. Through the comprehensive evaluation of the model indexes by the TOPSIS method with entropy value and t-test of model pairing， the fitting effect of these four age group dummy variable models （M20， M22， M25， M27） was better， which is suitable for the prediction of bark thickness of S. superba plantation in Fujian.

Keywords： bark thickness; dummy variable; age groups; site quality; Schima superba

DOI： 10.3969/j.issn.1000-2561.2020.12.027

树皮不仅对树木自身成长起到保护作用，还可成为生产及药用材料等重要生物资源，若能合理对树皮进行基质化利用还可进一步提高其经济价值[1]。树皮厚度是指树干某处带皮直径（diameter outside bark， DOB）与该处去皮直径（diameter inside bark， DIB）差的一半[2]。树皮大约占到整个木材体积的12%～20%，不同树种的树皮总量变化较大[3]，同一树种在不同位置的树皮厚度也不同[4]。因此，树皮厚度合理测算将直接影响到树皮蓄积量准确估算。

国外学者对树皮厚度的研究，早期主要集中在胸径处，如Grosenhbaugh[5]建立比率方程式，通过预测胸径高度处去皮直径来测算树皮厚度。Wingerd等[6]利用多元回归分析方法，构建美国阿巴拉契亚山脉7个树种的去皮直径回归模型，Maguire等[7]则利用削度方程来测算树干不同高度处去皮直径，这些方法都较大提高了树皮厚度计算的准确度。国内研究者通过加入胸径、树高、任意树高、年龄等解释变量，对白桦、落叶松、西南桦、湿地松、厚朴、云杉等树种的树皮厚度进行模拟，发现不同变量对树皮厚度影响也不同[8-13]。但目前研究还较少考虑到立地质量、龄组等因子对树皮厚度影响。贾炜玮等[14]建立含龄组、区域哑变量碳储量预测模型系统，结果显示含有哑变量的模型精度较高。朱光玉等[15]则构建含林分和立地类型哑变量的林分断面积生长模型，通过比较发现哑变量模型精度均优于基础模型。所以笔者将通过使用哑变量的方法将这些因子纳入到研究中。

木荷（Schima superba），也称为荷树、荷木，属于山茶科（Theaceae）常绿乔木[16-17]。木荷不仅是珍贵高级用材树种[18]，也是重要的防火树种[19]。近年来，福建和广东两省种植了大量的木荷人工林[20]，已有学者在木荷培育[21]、病虫害防治[22-23]等方面有研究，但还较少对木荷树皮厚度进行探讨。为此，本研究以福建省南平市木荷人工林为研究对象，通过对样木及样地调查获取模型所需因子，分别用胸高处树皮厚度（bark thickness at breast height， BBT）、任意高度处树皮厚度（bark thickness， BT）、相对树皮厚度（relative bark thickness， RBT）和去皮直径（diameter inside bark， DIB）4类基础模型[10]，拟合筛选出各自最优模型后，加入包含龄组和立地质量的哑变量，通过含熵值的TOPSIS综合分析法选出相容且统一的木荷树皮厚度模型，为计算木荷材积、出材量及生产经营措施制订提供参考依据，对木荷人工林经济价值精准评估具有重要现实意义。

1  材料与方法

1.1  材料

1.1.1  研究区概况  研究区为福建省南平市，位于福建北部，武夷山脉北段东南侧，全市面积约2.63万km2。南平地处北纬26°15～28°19，东经117°00～119°17之间，海拔50～2158 m。该地区属于中亚热带季风气候，年均气温19～22 ℃，无霜期300 d以上，年均日照时数1670～1980 h，年均降水量1400～1800 mm，相对湿度为79%。

1.1.2  数据获取  在福建省南平市木荷人工林中设置213块固定样地（含标准地），根据福建省森林资源清查中相关技术标准和实际林业生产情况，考虑不同海拔、地貌类型、林分年龄、胸径、树高、立地质量、郁闭度等因子，采取随机抽样的方式，选取189株样木，其中幼龄林19株、中龄林67株、近熟林49株、成熟林45株和过熟林9株;这些样木所在样地的立地质量等级分别为：（Ⅰ）肥沃15株，（Ⅱ）较肥沃68株，（Ⅲ）中等肥沃88株，（Ⅳ）贫瘠18株。根据国家气象科学数据中心提供的数据显示（http：//www.nmic. cn/），样木所在区域1981—2010年的极端气温、年均降水量、年均气温及年均湿度基本趋于一致，具体解析木所在标准地基本信息见表1。对样木伐倒后用精度为0.01 m的皮尺测量其全树高，打枝后以2 m为区分段，在树干的0.3、1.0、1.3、3.0、5.0 m等位置处用精度为0.01 cm的钢尺分别测量出带皮直径（DOB）和去皮直径（DIB）。

使用R语言简单随机抽样功能，从这189株样木中随机抽取两组数据，其中用于模型构建133株，其余部分用于模型检验（56株），进行本次研究模型的变量有：胸径（diameter at breast height， DBH）、树高（tree height， H）、林分年龄（stand age， T）、胸高处树皮厚度（bark thickness at breast height， BBT）、相对树皮厚度（relative bark thickness， RBT）、任意树高处树皮率（该处去皮直径与带皮直径之比，DIB/DOB）、郁闭度（canopy closure， CC），样木基本信息见表2。

1.2  方法

1.2.1  基础模型确立  本研究通过收集他人建立的14个模型，加上自建5个模型，共19个模型作为研究木荷树皮厚度基础模型（表3），探讨木荷最优树皮厚度模型。其中，树皮模型分类方法是根据唐诚等[10]的研究结果。

1.2.2  模型检验和评价  使用最小二乘法对以上19个基础模型进行拟合，计算各参数值，并进行显著性检验。经过筛选，对于满足所有参数均具有显著性的基礎模型，选取数学统计量指标：决定系数（determination coefficient， R2）、均方根误差（root mean square error， Erms）、和方差（the sum of squares due to error， Ess）、赤池信息准则（akaike information criterion， AIC）、贝叶斯准则（bayesian information criterion， BIC）等5个统计指标来评估基础模型拟合效果。判别标准为：R2越大，精度越高;Erms、Ess、AIC、BIC越小，精度越高;模型拟合效果也就越好。具体计算公式如下：

式中：yi为第i个实测值，i为第i个预测值，i为实测值的平均值，n为用于拟合模型的实测值数目，q为模型参数数目，为对数似然函数的最大值。

1.2.3  模型优选方法  TOPSIS法又可称为优劣距离法，1981年由Hwang和Yoon等[24]提出，它是一种有效的对多目标决策方法。但TOPSIS法在按欧氏距离进行优劣排序时，常出现接近于理想解，同时也可能出现贴近于负理想解的结果[25-26]。因此，需要用含熵值的TOPSIS法来弥补不足。熵值法是用来判断评价指标离散程度，指标离散程度越大则对综合评价影响也就越大[27-28]。具体计算步骤如下：

（1）根据模型求解结果构建m个评价对象、n个模型评价指标原始矩阵R。

（2）将原始矩阵进行归一化处理，得到无量纲化矩阵B，指标范围在0～1之间，值越大越优，最优值等于1，最差值为0。

xqmin为指标最小值，xqmax为指标最大值，rpq为指标归一化值，fpq为一致性指标，Hq为指标信息熵，ωq即为第q个指标的熵值，其中ωq不能大于1，lnn必大于0。

1.2.4  哑变量设置  哑变量，又称指示变量或虚拟变量，是对分类变量或定性因子进行处理的一种常用方法[29]。哑变量一般将定性因子转为定量因子，通常取值为0或1。在林业相关回归分析和数量化方法上经常引入哑变量，这不仅可以减少工作量，还可以提高計算准确度和模型精度[30-32]。本研究将通过优选基础模型，在不同类型树皮厚度模型中引入龄组和立地质量2个哑变量，用定性代码来表示，以构建具有相容性的模型[33]。

以龄组（i）做为哑变量为例，不同龄组使用定量变量Ji，表示如下：

当样木龄组为幼龄林时，取J1=1，J2=0，J3=0，J4=0，J5=0;当样木龄组为中龄林时，取J2=1，J1=0，J3=0，J4=0，J5=0;依此类推。同理，使用Ks表示立地质量，s=1，2，3，4，当样木来自立地质量为Ⅰ肥沃时，取K1=1，K2=0，K3=0，K4=0;当样木来自立地质量为Ⅱ较肥沃时，取K2=1，K1=0，K3=0，K4=0;依此类推。将Ji和Ks加入到各个模型参数中，以达到预测不同龄组和立地质量所包含的林分因子的目的。

1.3  数据处理

采用Excel 2010、SPSS 21.0、R软件（ver 3.6.1）进行各项数据处理分析。

2  结果与分析

2.1  基础模型拟合参数检验

将189株样木进行解析后进行模型拟合、参数求解并进行显著性检验后发现（表4）：在胸高处树皮厚度模型M2、M4、M5、M6中，均有1～3个参数预估值在5%水平与零差异不显著，且4个模型中都含有T、H等自变量，而最优模型M1中只含有胸径（DBH）自变量，说明胸高处树皮厚度与胸径更明显相关;在任意高度处树皮厚度模型M8、M9中，各有1个参数预估值在5%水平与零差异不显著，分别对应自变量为RH、CC、T;在相对树皮厚度模型M14、M15中，有1～2个参数预估值在5%水平与零差异不显著，分别对应自变量为RH、CC;而在去皮直径模型中未发现有参数不显著因子。因此这8个模型不再进一步参与分析，将剩下参数均显著的11个模型做进一步评价。

2.2  基础模型评价指标及优选

基于样木数据，对各基础模型进行拟合，结果显示：胸高处树皮厚度模型M1决定系数（R2 = 0.6981）在同类模型中最高;任意高度处树皮厚度模型M10决定系数（R2 = 0.5540）大于模型M7、M8;相对树皮厚度模型M12决定系数（R2 = 0.5056）大于模型M13;同时模型M1、M10、M12的Erms、Ess、AIC、BIC均为各类模型中最小（表5）。因此，模型M1、M10、M12分别为胸高处、任意高度处、相对树皮厚度最优基础模型。

由于去皮直径基础模型（M16、M17、M18、M19）的各评价指标互有优劣，因此，在对模型优选时，先对各评价指标进行无量纲化，再求解出各无量纲化值;然后利用无量纲化值计算各指标的熵值，其中Ess熵值最高为0.213102，R2熵值最低为0.161043;最后综合计算得出模型M16最优解总距离是0.000306，距离最小，为去皮直径最优基础模型（表6）。

2.3  含哑变量模型拟合参数检验和评价指标

为进一步研究模型，将龄组和立地质量哑变量分别同时加在基础模型的不同参数上，运用R软件进行分析与比较。结果显示，含龄组和立地质量哑变量模型的精度均要优于同时含有龄组和立地质量的哑变量模型。因此选出所有参数预估值与零差异显著模型进行评价。

由表7可知，在胸高处树皮厚度哑变量模型M20、M21的参数a中，分别加入龄组和立地质量哑变量后，R2从0.6981提高到0.9769和0.9758，但在Erms、Ess、AIC、BIC等指标上有略微增大，需进一步评价;在任意高度处树皮厚度哑变量模型M22、M23的参数a中加入龄组和立地质量哑变量后，其中模型M22的R2为0.9214，高于模型M23、M10，并且模型M22的Erms、Ess、AIC、BIC均小于同类其他模型，因此模型M22可直接进行模型t检验;在相对树皮厚度哑变量模型M24的参数a、b、c中加入龄组哑变量，以及在模型M25和M26的参数a中分别加入龄组和立地质量哑变量后，其中模型M25的R2最高为0.9111，但在Erms、Ess、AIC、BIC上3个哑变量模型和基础模型互有高低，需要进一步筛选;去皮直径哑变量模型M27、M28的参数b中，分别加入龄组和立地质量哑变量后，R2均为0.9997，但在Erms、Ess、AIC、BIC等指标上无法直接选出最优，需进一步分析。

Wingerd D E， Wiant H V. Variables for predicting inside-bark diameters of upper stems of appalachian hardwoods[J]. Journal of Forestry-Washington， 1982， 80（12）： 791-792.

Maguire D A， Hann D W. Bark thickness and bark volume in southwestern oregon douglas-fir[J]. Western Journal of Applied Forestry， 1990， 5（1）： 5-8.

王曉明. 兴安落叶松、白桦树皮厚度与胸径关系研究[J]. 林业科技情报， 2017， 49（4）： 45-47.

张如义， 邓小东， 杨小建， 等. 树高、胸径、树皮厚对厚朴产量的影响[J]. 经济林研究， 2017， 35（3）： 239-243.

唐  诚， 王春胜， 庞圣江， 等. 西南桦人工林树皮厚度模型模拟[J]. 林业科学， 2017， 53（7）： 85-93.

严铭海， 方静仪， 李睿宇， 等. 福州市湿地松人工林树皮厚度模型研究[J]. 西南林业大学学报（自然科学）， 2018， 38（5）： 124-131.

王晓林， 蔡可旺， 姜立春. 落叶松树皮厚度变化规律的研究[J]. 森林工程， 2011， 27（2）： 8-11.

郭孝玉， 孙玉军， 马  炜， 等. 适于FVS的长白落叶松树皮因子[J]. 东北林业大学学报， 2011， 39（10）： 28-31.

贾炜玮， 孙赫明， 李凤日. 包含哑变量的黑龙江省落叶松人工林碳储量预测模型系统[J]. 应用生态学报， 2019， 30（3）： 814-822.

朱光玉， 胡  松， 符利勇. 基于哑变量的湖南栎类天然林林分断面积生长模型[J]. 南京林业大学学报（自然科学版）， 2018， 42（2）： 155-162.

陈  林， 董安强， 王发国， 等. 广东南岭国家级自然保护区疏齿木荷+福建柏群落结构与物种多样性研究[J]. 热带亚热带植物学报， 2010， 18（1）： 59-67.

刘发林. 南方集体林区木荷次生林生长规律及经营技术研究[D]. 长沙： 中南林业科技大学， 2013.

阮传成， 李振问， 陈诚和， 等. 木荷生物工程防火机理及应用[M]. 成都： 电子科技大学出版社， 1995.

林  艳， 何紫迪， 毛积鹏， 等. 木荷基因组SSR位点开发及初步分析[J]. 热带作物学报， 2018， 39（9）： 1766-1771.

王秀花， 马丽珍， 马雪红， 等. 木荷人工林生长和木材基本密度[J]. 林业科学， 2011， 47（7）： 138-144.

黄  宇. 木荷无性系种子园生长与结实性状的初步研究[J]. 热带作物学报， 2017， 38（2）： 213-217.

邓志勇. 木荷提取物对小菜蛾和菜粉蝶生长发育的抑制效果[J]. 江苏农业科学， 2013， 41（4）： 129-130.

黄光强. 用木荷树皮粉末作为非杀伤性驱散武器填充剂的设想[J]. 森林公安， 2002（2）： 20.

Hwang C L， Yoon K. Multiple attribute decision making and applications[M]. New York： Springer-Verlag， 1981.

徐小万， 雷建军， 李  颖， 等. 现蕾期辣椒耐高温多湿性CA-TOPSIS综合评定[J]. 热带作物学报， 2013， 34（9）： 1747-1751.

Suthar H A， Gadit J J. Multiobjective optimization of 2dof controller using evolutionary and swarm intelligence enhanced with TOPSIS[J]. Heliyon， 2019， 5（4）： 321-332.

张一达， 刘学录， 范亚红， 等. 基于改进TOPSIS法的兰州市土地利用多功能性评价[J]. 干旱区地理， 2019， 42（2）： 444-451.

张浩渺. 基于熵权改进的TOPSIS模型在苏州河水质综合评价中的应用[J]. 西北水电， 2019（3）： 12-15.

曾伟生， 唐守正， 夏忠胜， 等. 利用线性混合模型和哑变量模型方法建立贵州省通用性生物量方程[J]. 林业科学研究， 2011， 24（3）： 285-291.

金  龙， 吴志祥， 杨  川， 等. 不同林龄橡胶凋落物叶分解特性与有机碳动态研究[J]. 热带作物学报， 2015， 36（4）： 698-705.

高东启， 邓华锋， 王海宾， 等. 基于哑变量的蒙古栎林分生长模型[J]. 东北林业大学学报， 2014， 42（1）： 61-64.

郑冬梅， 曾伟生. 用哑变量方法构建东北落叶松和栎类分段地上生物量模型[J]. 北京林业大学学报， 2013， 35（6）： 23-27.

唐守正， 郎奎建， 李海奎. 统计和生物数学模型计算： ForStat教程[M]. 北京： 科学出版社， 2009.

陈东来， 秦淑英. 树皮厚度、树皮材积与直径和树高相关关系的研究[J]. 河北林学院学报， 1994（3）： 248-250.

于世川， 张文辉， 李  罡， 等. 黄龙山林区不同郁闭度对辽东栎种群结构的影响[J]. 生态学报， 2017， 37（5）： 1537-1548.

Li R， Weiskittel A R . Estimating and predicting bark thickness for seven conifer species in the Acadian Region of North America using a mixed-effects modeling approach： comparison of model forms and subsampling strategies[J]. European Journal of Forest Research， 2011， 130（2）： 219-233.

池上评， 陈金章， 江传阳， 等. 基于间隔期的福建柏人工林动态生长模型及应用[J]. 福建林学院学报， 2014， 34（4）： 304-308.