湿地松(Pinus elliottii Engelm)原分布于古巴、美国东南部等区域^[1]，因其适应性强，干型通直，木材质量好^[2-3]，抗侵染能力优于马尾松^[4-5]，现已成为我国南方低丘陵地区的重要造林树种之一^[6-8]。湿地松的树皮用途很广，有研究表明其树皮提取物具有一定的药用价值，同时也是制作树脂胶的原料^[9]。树皮厚度可用于计算出材率、木材材积和树皮材积^[10]。国内外许多学者对树皮厚度进行了研究，涉及的树种有西南桦、银杉、云杉、白桦和厚朴等^[11-15]，有研究发现不同树种的树皮生长规律存在差异^[13]，通常针叶树的树皮厚度比阔叶树的大^[10]。对于树皮厚度模型的研究，拟合的因子主要有胸径、树高和年龄^[12]，但较少有考虑立地质量对树皮厚度的影响。如果考虑立地质量对树皮厚度的影响，则需对生长在不同立地质量等级的林木分别建立树皮厚度模型，在实际工作中这可能会耗费更多人力物力，且可能会出现模型不兼容的问题，因此采用哑变量模型可为兼顾模型拟合效果和成本提供途径^[16]。

已有研究将立地质量等级、竞争状态、经营管理措施等定性因子作为哑变量引入到模型中，使模型的拟合效果得到提升。朱光玉等^[17]在研究栎类天然林林分断面积生长模型时引入了立地质量等级为哑变量，华伟平等^[18]在黄山松地位级指数模型的研制中引入立地质量等级为哑变量，曹梦等^[19]在单木胸径和树高生长模型研究中引入不同竞争状态为哑变量，ZENG^[20]在生物量的研究中引入森林起源为哑变量，王金池等^[16]在云南松林分蓄积量生长模型的研究中引入间伐与未间伐指标为哑变量。因此本研究通过引入立地质量等级作为哑变量，构建湿地松树皮厚度哑变量模型，以达到提升模型精度的目的。

本文通过参考前人研究所得的树皮厚度模型，利用湿地松样地和样木数据，得到了湿地松带皮直径、去皮直径和相关的树皮厚度数据，利用含熵权值的TOPSIS综合分析法选出最优的树皮厚度模型，在此基础上，引入立地质量作为哑变量建模，构建了基于哑变量的湿地松树皮厚度模型，为测算湿地松木材材积、出材率以及树皮蓄积量提供依据。

1 研究区概况

宁德市地处福建省东北部(东经118°32′~120°44′，北纬26°18′~27°04′)，土地面积13 400 km²，以丘陵山地为主兼沿海小平原, 属于亚热带海洋性季风气候，温暖湿润，降水充足。年平均气温13.4~20.2 ℃，年平均降雨量1 250~2 350 mm, 无霜期270.4 d，日照时间1 637.7 h。

2 研究方法 2.1 材料来源

在宁德市湿地松人工林中设置样地共86块，包括临时样地与固定样地，样地形状为方形与矩形，面积为0.04~0.67 hm²，设置于不同立地质量等级的林分中。结合伐区设计数据，从中选取涵盖不同年龄、胸径与立地质量的样木258株，所在样地基本信息见表 1。伐倒后用皮尺(精度0.01 m)测量树高与冠长，打枝后以2 m为区分段，用钢尺(精度为0.1 cm)于树干上0.3、1.0、1.3、3.0、5.0 m处依次测量带皮直径(outside bark diameter，dob)和去皮直径(inside bark diameter，dib)。随机抽取样木193株(约占75%)用于构建模型，剩余65株用于模型检验(约占25%)。样木数据的描述性统计结果见表 2，模型变量包括林分年龄(tree age，t)、树高(tree height，H)、胸径(diameter at breast height，DBH)、胸高处树皮厚度(bark thickness at breast height, BBT)、任意树高处树皮率(该处两倍树皮厚度与带皮直径的比值，Bark ratio at any height，z)、为任意树高处树皮厚度(bark thickness at any height，BT)、相对树皮厚度(任意树高处树皮厚度与胸径处树皮厚度的比值，relative bark thickness，RBT)、树冠长(crown length，CL)和树冠率(冠长与树高的比值，crown rate，CR)等。

2.2 哑变量定义

哑变量(虚拟变量)是定性变量，通常取值为0或1，这种方法称为定性因子(0, 1)化展开，即将变量δ(x，i)表示为:δ(x，i)=0或1，当x是第i等级时为1，否者为0, 因此称变量δ(x，i)为哑变量。本研究引入立地质量等级作为哑变量，根据福建省森林资源清查中的相关技术标准和实际林业生产情况将立地质量分为Ⅰ肥沃、Ⅱ较肥沃、Ⅲ中等肥沃、Ⅳ贫瘠4个等级。将立地质量等级用定性代码0或1表示，第i种立地质量等级表示为S_i，定性数据S_i转化为(0，1)形式:S_i为0或1，当x时第i等级时为1，否者为0。其中i=1、2、3、4；S₁、S₂、S₃、S₄分别是Ⅰ肥沃、Ⅱ较肥沃、Ⅲ中等肥沃、Ⅳ贫瘠的定性代码。

2.3 树皮厚度模型选择

通过查阅相关资料，参考唐诚等^[21]、严铭海等^[22]建立的17个模型，并用于本研究，模型表达式见表 3。其中模型3、模型4和模型5属于对数回归模型，需要计算校正因子进行修正，本文采用目前应用较多的校正因子(correction factor，CF)^[23]，其计算公式如下：

$ {C_{\rm{F}}} = {\rm{exp }}({s^2}/2) $

(1)

式中：C_F为校正因子；s²为回归方程样本方差。

2.4 模型评价指标

采用最小二乘法对17个模型进行拟合，对模型参数进行显著性分析，选取所有参数均达到显著的模型，采用赤池信息准则(Akaike information criterion，AIC)、贝叶斯准则(Bayesian information criterion，BIC)、偏差(deviation，B)、绝对偏差(absolute deviation，AB)和决定系数(decision coefficient，R²)5个指标做模型质量综合评价，计算公式如下：

$ B = \sum\limits_{i = 1}^n {\frac{{{y_i} - {{\hat y}_i}}}{n}} $

(2)

$ {A_{\rm{B}}} = \sum\limits_{i = 1}^n {\frac{{{y_i} - {{\hat y}_i}}}{n}} $

(3)

$ {R^2} = 1 - \frac{{\sum\limits_{i = 1}^n {{{({y_i} - {{\hat y}_i})}^2}} }}{{\sum\limits_{i = 1}^n {{{({y_i} - \bar y)}^2}} }} $

(4)

$ {B_{{\rm{IC}}}} = - 2{\rm{ln}}\ell + {\rm{ln}}np $

(5)

$ {A_{{\rm{IC}}}} = - 2{\rm{ln}}\ell + 2p $

(6)

式中：A_B为绝对偏差；B_IC为贝叶斯准则；A_IC为赤池信息准则；y_i为实测值；$\hat y$为预测值；$\bar y$为实测值的平均值；n为用于拟合模型的实测值数目；p为模型参数数目；ln $\ell $为对数似然函数的最大值。

2.5 模型优选方法

本研究采用含熵权值的TOPSIS法，在计算模型评价指标权重时运用熵值法取代一般的主观权重法，避免了人为主观性的影响^[24-26]，其计算步骤如下。

(1) 根据模型指标求解结果构建m个评价对象、n个评价指标的判断矩阵R，公式如下：

$ ~\mathrm{R}={{({{x}_{ij}})}_{m\times n}}(i=1, 2, \ldots , m;j=1, 2, \ldots , n)$

(7)

(2) 将判断矩阵归一化处理，得到无量纲化后的矩阵B，求解出的指标值范围为0~1，其值越大代表指标越优，最优解的值为1，最差解值为0。

高优指标(效益型指标)公式如下：

$ {{r}_{ij}}=\left\{ \begin{align} & \frac{{{x}_{ij}}-{{x}_{j, \min }}}{{{x}_{j, \text{max}}}-{{x}_{j, \text{min}}}}({{r}_{j, \text{max}}}\ne {{r}_{j, \text{min}}}) \\ & 1({{x}_{j, \text{max}}}={{x}_{j, \text{min}}}) \\ \end{align} \right.\text{ } $

(8)

低优指标(成本型指标)公式如下：

$ {{r}_{ij}}=\left\{ \begin{align} & \frac{{{x}_{j, \text{max}}}-{{x}_{ij}}}{{{x}_{j, \text{max}}}-{{x}_{j, \text{min}}}}({{r}_{j, \text{max}}}\ne {{r}_{j, \text{min}}}) \\ & 1({{x}_{j, \text{max}}}={{x}_{j, \text{min}}}) \\ \end{align} \right.\text{ } $

(9)

(3) 计算熵权值W_j

$ {{H}_{j}}=-\frac{1}{\text{ln}~n}(\sum\limits_{i=1}^{n}{{{f}_{ij}}\text{ln}~{{f}_{ij}}}~)(i=1, 2, \ldots , m;j=1, 2, \ldots , n) $

(10)

$ ~{{f}_{ij}}=\frac{1+{{r}_{ij}}}{\sum\limits_{j=1}^{n}{({{r}_{ij}}+1)}} $

(11)

$ {{W}_{j}}=\frac{1-{{H}_{j}}}{n-\sum\limits_{j=1}^{n}{{{H}_{j}}}}~\left[ W={{({{W}_{j}})}_{1\times n}}, \sum\limits_{j=1}^{n}{{{W}_{j}}=1} \right] $

(12)

式中：x_ij为模型评价指标值；x_j，max为指标最大值；x_j，min为指标最小值；r_ij为评价指标归一化值；n为评价指标个数；f_ij为一致性指标；H_j为指标信息熵；W_j为熵权值。

3 结果与分析 3.1 树皮厚度模型拟合结果

本文利用258株样木数据求解出17个树皮厚度模型并对模型参数进行显著性检验(表 4)。结果表明，模型2、模型5、模型7、模型11均有一个参数不显著，分别对应的自变量为H、ln H、t、H_R²，将这4个模型予以剔除。其中H和ln H是胸高处树皮厚度模型的因子，t是任意高度处树皮厚度模型的因子，H_R²是相对树皮厚度模型的因子，去皮直径模型中没有出现参数未显著的因子。经模型参数显著性检验后，筛选出14个参数均显著的树皮厚度模型，这些模型进入下一步的模型质量评价。

3.2 模型质量评价

模型拟合效果与预估精度直接相关，选用的5个评价指标从不同方面体现模型的优度，分别对4类模型的评价指标进行综合分析，模型指标计算结果及含熵权值TOPSIS法结果见表 5，各评价指标熵权值见表 6。

在模型评价指标无量纲化后，指标取值范围为0~1，其值越接近1代表指标优度越高，其值越接近0代表指标优度越低。模型3评价指标的B值在胸高处树皮厚度模型中最小为0.000 002，该指标无量纲化后取值为1，代表该模型的评价指标B优度在胸高处树皮厚度模型中最高。模型6评价指标的A_IC值在任意高度处树皮厚度模型中最大为6 121.3，该指标无量纲化后取值为0，代表该模型的评价指标A_IC优度在任意高度处树皮厚度模型中最低。

综合本文所选的各评价指标无量纲化值及各评价指标的熵权值得到的结果表明，胸高处树皮断面积模型相对最优解距离分别为0.194 3>0.034 3>0.000 9，表明模型质量模型3>模型4>模型1，评价指标中A_B熵权值最高为0.211 0，B熵权值最低为0.188 9。任意高度处树皮厚度模型相对最优解距离分别为0.115 3>0.086 5>0.000 4>0.000 0表明模型质量模型9>模型10>模型8>模型6，评价指标中AB熵权值最高为0.206 3，R²熵权值最低为0.194 2。相对树皮厚度模型相对最优解距离分别0.200 4>0.042 2>0.000 4表明模型质量模型12>模型13>模型14, 评价指标中AB熵权值最高为0.215 1，B熵权值最低为0.190 6。去皮直径模型相对最优解距离分别0.001 5>0.000 4>0.000 0表明模型质量模型16>模型15>模型17, 评价指标中AB熵权值最高为0.205 0，R²熵权值最低为0.185 0。由此，选出模型3、模型9、模型12、模型16这4个相对最优解距离最小的最优模型。

3.3 哑变量模型及其拟合结果

在这4个最优模型引入哑变量，经过尝试发现，在不同参数中加入哑变量时模型的决定系数有差异，模型3、模型9、模型12、模型16分别在参数α、α、β₂、α中加入哑变量后所得模型不存在未显著参数且拟合效果最佳，因此，哑变量模型的形式见表 7。其中b_i为哑变量S_i的参数(i=1，2，3，4)，当立地质量为Ⅰ肥沃时S₁=1、S₂=0、S₃=0、S₄=0，Ⅱ较肥沃时S₁=0、S₂=1、S₃=0、S₄=0，Ⅲ中等肥沃时S₁=0、S₂=0、S₃=1、S₄=0，Ⅳ贫瘠时S₁=0、S₂=0、S₃=0、S₄=1。以建模数据拟合以上4个哑变量模型，其参数求解值与各评价指标值见表 8和表 9。

由表 9可知，4类哑变量模型的参数均显著，4类哑变量模型较对应的最优基础模型决定系数R²均得到提升，其中模型18的R²为0.592，模型19的R²为0.862，模型20的R²为0.868，模型21的R²为0.997。哑变量模型较最优基础模型A_IC、B_IC等低优指标均有降低，其中模型18为341.3、356.7；模型19为5 045.7、5 102.2；模型20为919.7、955.4；模型21为13 011.1、13 053.3，说明所建的哑变量树皮厚度模型的拟合效果在筛选出的最优模型的基础上得到了提升。

3.4 模型拟合效果

根据用于检验的65株样木实测数据，来进一步检验哑变量模型对生长于4种立地质量等级的湿地松树皮厚度预估效果。采取配对t检验的方法检验模型拟合效果，假设模型预估值与实测值之间存在显著差异。若P＜0.05，则假设成立，表示模型的拟合效果不好；若P>0.05，则假设不成立，表示模型拟合效果优良(表 10)。

结果表明，4类湿地松树皮厚度哑变量模型的P值均大于0.05，其值分别为0.811、0.955、0.805、1.000。配对t检验结果表明本文所建立的4类湿地松树皮厚度哑变量模型的预估值与实测值之间无显著差异，表明所建模型可以做为计算生长于Ⅰ肥沃、Ⅱ较肥沃、Ⅲ中等肥沃、Ⅳ贫瘠4种不同立地质量下湿地松树皮厚度的参考。

4 讨论与结论

本研究对拟合出的模型进行显著性检验，剔除含不显著参数模型后运用5个评价指标进行综合评价。依据含熵权值TOPSIS法评价结果筛选出最优模型3、模型9、模型12、模型16，在这4个模型的基础上引入哑变量，4个哑变量模型较基础模型拟合效果均得到提升。哑变量模型经配对t检验后，其P值分别为0.811、0.955、0.805、1.000，检验结果表明4类哑变量模型对检验样本的预估值与其实测值不存在明显差异，可以做为计算不同立地质量下湿地松树皮厚度的参考。对4类树皮厚度哑变量模型的S_i参数值进行分析，发现模型18的哑变量的S_i参数值b₄ < b₃ < b₂ < b₁, 说明生长于不同立地质量的湿地松胸高处树皮厚度由大到小的排序为Ⅰ肥沃>Ⅱ较肥沃>Ⅲ中等肥沃>Ⅳ贫瘠；模型19的哑变量的S_i参数值b₄ < b₃ < b₂ < b₁, 说明生长于不同立地质量的湿地松任意高度处树皮厚度由大到小的排序为Ⅰ肥沃>Ⅱ较肥沃>Ⅲ中等肥沃>Ⅳ贫瘠；模型20的哑变量的S_i参数值b₂ < b₄ < b₁ < b₃, 说明生长于不同立地质量的湿地松相对树皮厚度由小到大的排序为Ⅲ中等肥沃>Ⅰ肥沃>Ⅳ贫瘠>Ⅱ较肥沃；模型21的哑变量的S_i参数值b₄ < b₃ < b₂ < b₁, 说明生长于不同立地质量的湿地松去皮直径由大到小的排序为Ⅰ肥沃>Ⅱ较肥沃>Ⅲ中等肥沃>Ⅳ贫瘠。其中模型18、模型19、模型21均是在常数项加入哑变量的拟合效果最佳，依此可见立地质量对湿地松胸高处树皮厚度、任意树高处树皮厚度和去皮直径可能有比较直接的影响而且是加法关系。模型20在常数项加入哑变量后出现了参数不显著的情况，在指数项引入哑变量后模型的拟合效果最好。

基础模型选优采用含熵权值的TOPSIS法，选用的各项评价指标的熵权值由其本身的值决定，该方法具有较好的客观性^[24]，评价过程避免了人为主观的确定各个指标权重而影响综合评价的结果。本文中4种类型的模型各项评价指标熵权值均有不同，其中不同指标在同一类型模型中的重要性不一致且同一指标在不同类型模型中的重要性也存在差异，因此有必要采用改进的含熵权值的TOPSIS法对模型进行综合评价。

本研究考虑到立地质量对树皮厚度的影响，从理论角度，将立地质量等级作为哑变量加入到模型的构建中，能规避由于立地质量不同而出现模型预测有偏和不同立地质量单独建模不相容的情况。从实际工作的角度，在模型中引入哑变量，能避免重复建模问题且降低了调查的难度从而减小工作量。同一模型对不同参数引入哑变量得到的模型其决定系数存在不同，在代表树皮厚度的最大参数值中引入立地质量等级哑变量模型的预估效果最好，模型引入哑变量后决定系数增大，各项低优指标值均降低，模型的拟合效果得到了提升。模型的选用不仅要考虑准确性还需要实际操作的简易性，本研究根据唐诚等^[21]和严铭海等^[22]的研究选取了调查比较容易获得的因子。以胸径为主要调查因子加入了其他因子来提升模型的拟合效果，在拟合结果中出现了不显著的参数因子分别为H、ln H、H_R²、t，其中H和ln H因子在胸高处树皮厚度模型中出现了不显著的参数，t因子在任意高度树皮厚度模型中出现不显著参数, H_R²在相对树皮厚度模型中出现了不显著参数。在树高相关的因子中出现不显著的参数与前人^[21]的研究相符合。本文通过样木和样地的数据初步构建了4类含立地质量等级哑变量的树皮厚度模型，为测算湿地松出材率及树皮蓄积量提供依据，但还没有考虑到气候因子对树皮厚度的影响，往后的研究可以尝试将气候因子引入模型中来提高精度。