树冠大小与树木的光合能力密切相关^[1]，直接影响着树木生活力和生产力，代表了过去一段时间内树木的累积竞争效应^[2]，其中冠幅(crown width，CW) 作为衡量树冠大小的主要指标之一，通过预测冠幅生长可用来计算树冠的表面积和体积等^[3]，是林分可视化的重要参数^[4]。此外，冠幅也可以用来表示单一林层中树木之间的竞争程度^[1]，然而在异龄复层的天然林中，林层结构复杂，上层树木树冠会遮挡其冠下树木对太阳直射光的吸收^[5]，此时，树冠活枝下高(height to crown base，HCB，指树冠的第一活枝到地面的高度) 作为反映树冠所处垂直结构动态变化的重要指标，直接反映了树冠在垂直方向的生长和衰退，此外，枝下高可以用来计算冠长、冠长率等树冠大小^[6]指标，其作为单木预测因子也常用于预测冠幅模型^[7]、冠型模型^[8]和生物量模型^[9]等。可见，树木枝下高和冠幅有巨大应用价值，但是测量每株树木的枝下高和冠幅会浪费大量的人力物力，且在高密度的复层林中，枝下高和冠幅的测量精度差，因此研究树木枝下高^[10]和冠幅^{[4, 7]}的预测模型具有重要的实际意义。

树木的光合作用影响着树冠冠幅大小^{[3, 4, 8]}，随着树木生长，单木下部枝条受到的自身和相邻木的遮阴逐渐增加，当光照强度低于光合作用光补偿点时，树木进行自然整枝^[11]，使得枝下高增高，冠幅减少。不同树种的光补偿点和耐荫性不同使得其树冠衰退速率不同^[12]。此外，对于单木大小相近的树木，竞争则会促进林木的高生长，也加剧了相邻木树冠间的遮阴，使得林木有更高的枝下高^[11]和更小的冠幅^[4]；贺梦莹等^[7]在长白落叶松-水曲柳混交林的冠幅模型中发现大于对象木胸高断面积之和、单木胸径与林分优势胸径之比等与距离无关的竞争因子与冠幅显著相关。同样的，立地因子中坡向和坡度差异使得林地的日光照时长和角度不同，进而影响单木树冠之间的遮阴效应，国内外学者在对枝下高^{[10, 12]}和冠幅^[13]的模型研究中，发现添加了立地因子的模型拟合效果有了显著提高^[9]。

为了解单木的树冠生长规律，不少学者展开了相关模型研究。已有研究多集中在对树冠的单一指标进行回归分析，如：雷相东等^[4]采用线性回归建立了落叶松云冷杉天然林内多树种的冠幅模型，卢军等^[14]建立了帽儿山天然次生林主要树种冠长率模型。也有学者以样地为随机效应，对枝下高^[6]或冠幅^{[7, 13]}构建混合效应模型以提高模型的预估精度。周泽宇等^[15]采用多种加权的参数估计方法构建了云冷杉更新幼树的冠幅可加性模型，验证了加权似不相关回归估计的参数估计法优于加权最小二乘法估计。李想等^[12]则应用联立方程组对人工樟子松的树高和枝下高模型进行系统估计，提高了模型的预估效果。然而在树冠指标中，冠幅反映了树木向周围空间的水平延伸长度，常用来表示树木的生长效率及其在林分中的竞争能力^[1]，而枝下高反映了林木营养空间的高度动态变化，树冠枝下高和冠幅具有高度内在相关性^{[4, 7-8]}，在天然异龄复层的混交林中，树冠活枝下高和冠幅分别从垂直和水平角度反映了单木树冠的竞争状态，通过对单木同时构建枝下高和冠幅模型，能够对天然林树冠大小和树冠衰退规律有更加完备的认识。

云冷杉林也称暗针叶林，是由云杉属(Picea) 和冷杉属(Abies) 树种组成的常绿、阴暗、潮湿的森林群落总称，我国云冷杉资源极为丰富，面积大、蓄积高，是东北地区重要的森林资源。本研究基于天然云冷杉针阔混交林内云杉(Picea jezoensis var. microsperma)、冷杉(Abies nephrolepis)、红松(Pinus koraiensis)和林内阔叶树种枝下高和冠幅数据，采用最小二乘法估计(least square method，OLS)、极大似然估计(maximum likelihood estimation，MLE) 和似不相关回归估计(seemingly unrelated regression estimation，SUR) 3种参数回归方法构建单木的枝下高和冠幅模型，对其拟合和检验精度进行对比，选出拟合和预测最优的模型来分析树种、单木大小、竞争和立地等因子对枝下高和冠幅模型的影响，为提高树冠生长模型预测精度和科学经营提供参考。

数据来源于吉林省汪清林业局金沟岭林场(43°17′~43°25′N，130°05′~130°20′E)，在云冷杉针阔混交林生态系统中选择有代表性的林地，面积92.5 hm²，海拔600~780 m，林地平均坡度15°，属温带大陆性季风气候，年平均气温3.9 ℃，年平均降水量600~700 mm，多集中在7月，主要的土壤种类为暗棕壤。研究区内主要树种有云杉、冷杉和红松，其他树种还有长白落叶松(Larix olgensis)、白桦(Betula platyphylla)、枫桦(Betula costata)、大青杨(Populus ussuriensisi)、蒙古栎(Quercus mongolica)、紫杉(Taxus cuspidata Zucc.var. chinensis)、春榆(Ulmus japonica)、色木槭(Acer pictum Thunb. ex Murray)、水曲柳(Fraxinus mandshurica)、黄檗(Phellodendron amurense) 和胡桃楸(Juglans mandshurica) 等。

在研究区内机械设置98块长期固定样地，每块样地0.04 hm² (20 m×20 m)，样地间距90 m，在每个样地中心埋设一个水泥标桩，注明样地号，并将样地内的每棵树钉牌，用于编号，以便于复查。2019年7月，对98块样地进行调查，记录样地海拔、坡度和坡向(其中北记为0°，西北记为45°，西为90°，以此类推) 等立地因子，以及林分内胸径(diameter at breast height，DBH) ≥5.0 cm的树木胸径、树高(tree heigh，H)、枝下高和冠幅等因子。样地的林分特征和立地因子统计如表 1所示，其中在长期固定样地布设过程中，严格按照等距离的机械布置，该研究区经常受台风影响，林区内高大树木的风折风倒较为常见，反映在固定样地中则为极少部分样地受台风风折影响，样地损坏严重，样地内胸径≥5 cm的活立木较少，为了更为真实地反映研究区内的林分现状，本次研究对于样木较少的样地给予保留。

样地中树种组成主要以云杉、冷杉和红松三大针叶树种为主，同时还分布有少量阔叶树。考虑到树种生长速率^[16]的不同，本研究将阔叶树种分为中生阔叶树种组(moderate growing hardwoods，下文用Hm表示，包含春榆、白桦和枫桦以及极少量的大青杨等) 和慢生阔叶树种组(slow growing hardwoods，下文用Hs表示，包含水曲柳、椴树、蒙古栎、色木槭和黄檗等)。共测得样木2 487株，其中云杉298株，冷杉646株，红松294株，中生阔叶281株，慢生阔叶968株。分树种(组) 将数据按照7∶3的比例随机划分为建模数据和检验数据，各树种基本单木测量指标统计见表 2。

本研究将变量分为4组：单木大小、林分因子、单木竞争和立地条件，各候选自变量及说明见表 3。对模型入选自变量采用方差膨胀因子(VIF) 来判断自变量间的多重共线性，仅保留VIF < 5且回归系数显著的自变量^{[4, 7]}。

胸径是预测冠幅的重要影响因子^[3]，而枝下高也是预测冠幅的重要因子^[4]，枝下高基础模型常见的有指数模型和Logistic模型，而冠幅常见的有生长型模型、多项式模型、Logistic模型和指数模型，在将拟合数据初步代入后，通过比较模型的修正决定系数和均方根误差等评价指标，最终采用Logistic模型^{[10, 17]}和指数模型^[1]来分别预测枝下高和冠幅的生长规律：

式中：a₁~a₄、b₁、c₁~c₂为参数；b₀、c₀为常数项；V_i为参数组，其中V₁为单木大小组，V₂为林分因子组，V₃为单木竞争组，V₄为立地条件组。2个模型的同一次观测模型误差ε₁和ε₂的各分量间是相关的，即cov (ε₁，ε₂) 是非对角矩阵。

本研究最小二乘法估计和极大似然估计是对枝下高和冠幅模型分别进行拟合，而似不相关回归估计是对同一对象木的枝下高和冠幅2个模型共同进行拟合，使得每个树种(组) 的2个方程拥有同一套参数。数据分析主要采用R 4.0.3软件，使用systemfit包构建最小二乘法估计和似不相关回归估计模型，nlme包构建极大似然估计模型^[18]。

本研究基于模型的基础形式，其拟合和检验的结果采用以下指标评价^[4]：修正决定系数(R_a²)、均方根误差(E_RMS) 和平均绝对误差(E_MA)。最优模型选择根据R_a²最大，E_RMS和E_MA最小的原则进行，其计算公式如下：

式中: E_RMS为均方根误差；y_i为实测值；${{{\hat y}_i}}$为预估值；y_i为实测的平均值；n为样本个数；p为模型参数个数；E_MA为平均绝对误差。

由表 4可得，3种参数回归估计方法的模型拟合R_a²、E_RMS和E_MA无明显差异；OLS模型拟合效果较SUR和MLE表现略好。而用检验数据对OLS、SUR和MLE 3种方法构建的枝下高和冠幅模型进行检验，发现采用SUR法的模型检验R_a²和E_RMS在各模型种中优于OLS和MLE；且所有模型检验的E_MA均较OLS和MLE降低了57%~58%。

通过3种参数估计模型的拟合和检验统计量比较，采用似不相关回归估计对枝下高-冠幅联立方程组模型的参数估计模型检验的E_MA较OLS和MLE降低了57%~58%，模型预测效果最好，因此以下分析均基于SUR法的参数估计结果。从表 5看出，所有树种的枝下高均与胸径呈正相关；由式(2) 得，枝下高与树高呈正相关，且枝下高随树高的增长趋势随着树高的增加而减缓；除红松外，所有树种均与单木相对胸径大小(D_BD) 呈显著负相关；除云杉外，其余树种枝下高均与样地株数(N) 显著负相关，红松枝下高与树冠竞争因子(C_F) 显著负相关，所有立地条件因子中，仅冷杉枝下高与反映样地光照时长的坡向余弦值(A_C，该值越大则样地受到的光照时长越少) 显著正相关。

所有树种冠幅与胸径(D) 和样地树冠竞争因子(C_F) 均显著正相关；云杉、冷杉和红松3种针叶树种的冠幅均与枝下高显著负相关，阔叶树种冠幅则与冠长(L_C) 显著正相关；云杉冠幅与样地断面积之和(B_A) 显著负相关；而冷杉、中生阔叶和慢生阔叶的冠幅则与样地株数(N) 呈负相关。

在表 5的参数估计中，N、D_g和C_F这3个林分因子为单木枝下高和冠幅模型的主要解释变量，因此，研究中假定一个样地环境，其N、D_g和C_F取所有样地平均值，用胸径^[12]或树高^[10]来模拟单木大小变化。图 1 (a)为在相同样地环境(N=29，D_g=18.5，C_F=140) 中云杉、冷杉、红松、中生阔叶和慢生阔叶5个树种(组) 的枝下高与胸径的关系曲线。图中，各个树种枝下高的增加速度先快后慢，针叶树种中，云杉的枝下高在同胸径大小的单木中最大，冷杉次之，红松的最低；两种阔叶树种的枝下高变化差异较小。同理，图 1(b)为5个树种(组) 的冠幅与胸径的关系曲线，与枝下高不同，单木冠幅的增长速率随着胸径的增加而增加；针叶树种中，同胸径大小的红松冠幅最大，其次是云杉，冷杉最小；同胸径大小的中生阔叶的冠幅远大于慢生阔叶。

	注：曲线横坐标范围为各树种胸径观测值区间。 Notes: The range of fitting curve is the observation interval of DBH of each tree species 图 1 5个树种(组) 的枝下高和冠幅与胸径的关系曲线 Fig. 1 Relationship between height to crown base and crown width with the diameter at breast height of 5 species
图选项

本研究应用SUR、OLS和MLE 3种方法对模型进行了参数估计，在对模型的拟合统计量R_a²和E_RMS比较中，OLS表现最好^[19]，这是由于OLS能够通过最小化误差的平方和寻找数据最佳函数进行匹配，这与周泽宇等^[15]在比较参数估计法差异研究中得到的结果一致，即SUR模型拟合R_a²和E_RMS比OLS略低，但相差不大。OLS作为回归模型的基础参数估计方法，在林业中的材积表绘制、温度模拟等方面应用较为广泛。根据SZATROWSKI^[20]和吴密霞^[21]的研究发现，固定效应的极大似然估计的模型效果接近其最小二乘估计，本研究中也得出了这一结论。最小二乘估计由于没有利用误差项方差－协方差矩阵的结构信息，其估计效率往往较低，在近年来的混合效应模型研究中，研究者主要采用极大似然估计的方法^{[6, 13]}，且在添加了模型随机效应后，预测效果得到显著提升^[3]，在特定林区内能发挥较好的预估效果，然而混合效应模型中随机效应参数受样地限制，不具备广泛区域适用性。而利用SUR法对单木枝下高和冠幅联立方程组模型的参数进行同时估计，不仅确保了2个方程参数估计的一致性，同时在对模型检验时，模型的E_MA较OLS和MLE显著降低^[19]，且就模型检验R_a²和E_RMS而言，SUR法同样在三者中表现最佳。在林业模型研究中，对于存在多个因变量的联立方程组模型，当变量间存在高度内相关性时，采用似不相关回归估计不仅考虑了方程组的同期偏差和同期相关性，还能够提供更加精准和稳定的参数估计值^{[12, 15]}。

本研究中枝下高的基础模型为Logistic形式的Wykoff模型^[17]，如图 1可得，所有树种的枝下高^[6]和冠幅^[13]均随着单木个体大小的增加而增加，且增长趋势逐步减缓，这与YANG et al^[10]在模拟白云杉、樟子松和白杨枝下高随树高的变化曲线中的结果相一致。树木在个体较小时，在林地中受到其他树木的遮阴较多^[7]，为了竞争到更多的光资源，优先进行高生长^[22]，随着树高的增加，在林分中竞争能力较弱的树木被自然稀疏后，林分密度下降，透光率增加^[23]，高大树冠衰退速率降低，因此，单木树冠衰退速率与个体大小呈负相关^[14]。此外，小树树冠下部枝条频繁地自我整枝，导致树木冠幅增长较慢^[8]。而当树木进入主林层后，树高受树种和环境限制接近极限高度后无法增加^[24]，自身在林内受到的竞争压力较小，枝下高将趋于稳定^{[12, 17]}，高大树木为了维持其竞争优势，将更多的营养位置用于树冠枝条的生长^[5]，进而截获更多的光资源，通过模拟枝下高和冠幅分别与胸径的关系曲线时也得出了相同的结论。在所有单木竞争指标中，树木胸径与样地平均胸径之比作为衡量树木在样地中相对大小的指标，是本研究中所有树种枝下高模型中较为一致的预测因子。单木在样地中的相对地位越低，其树冠受到更多相邻木的遮阴，树冠下层枝条整枝强度越高^[11]，导致枝下高越高^[25]。在图 1中，红松较同胸径大小的云冷杉拥有更低的枝下高和更大的冠幅，研究发现红松与反映样地内潜在立地质量的树冠竞争因子显著负相关^[26]，在立地质量较高的林地中，高大红松在生长过程中能够从环境中获取更多的光照、水和养分，进而维持其树冠下部粗大枝条的生长^[27]，使其树冠衰退速率变低。而其余树种在树冠竞争因子大的样地中，能获取更多营养物质，有利于冠幅的生长^[13]，但也会受到更多上层树木的树冠遮阴，导致其枝下高与树冠竞争因子关系不显著。

云杉和冷杉为该地区常见的伴生树种^[28]，但由于云杉比冷杉的耐荫性相对较弱^[29]，树冠下部枝条整枝强度更高，使得林分中胸径相近的云杉比冷杉的枝下高更高(图 1)。在立地条件较好的且以云冷杉为主的样地中，如表 2所示，尽管冷杉在顶层林未能得到很好的生长空间，但是在林下个体较小的更新幼树中，通常多为冷杉^[28]。由于冷杉和阔叶树在研究区内分布密集且个体较小，因此林分密度能够很好地反映其相互间的竞争强度，在林分密度更大的样地中林木会经历更多的竞争，竞争的加剧使得树冠衰退，枝下高增高^{[10, 27]}，冠幅减小^{[7, 13]}。在该森林类型中，云杉在林分的自然更新较少^[28]，云杉平均胸径比冷杉大10 cm，其株数不足冷杉的二分之一(表 2)，因此，对于低密度样地下的高大云杉，株数密度不能更好地反映其树木的竞争激烈程度，且样地中个体较大的云杉冠层已超过林下密集分布的小径阶冷杉^[29]，小个体的树木对高大树木的光竞争几乎可以忽略^[5]，导致云杉枝下高和冠幅均与所处样地的株数密度关系不显著。在云杉冠幅模型中，冠幅与林分断面积和显著负相关，林分断面积和相较于株数密度，更能够反映不同密度的林分中的云杉所受到的竞争压力。在添加样地坡向、坡位和海拔等立地因子进入树种枝下高模型中，发现冷杉枝下高与样地的每日光照时长负相关，树木在低密度群落中对光资源的竞争较小^[11]，而在高密度的中小径阶冷杉群落中，冷杉枝下高对光照差异的反应较为灵敏^[30]。由图 1可知，该森林类型中两组阔叶树种的枝下高差异较小，意味着两种树种组间耐荫性和树冠衰退速率差异较小，中生阔叶相较于同个体大小的慢生阔叶有着更大的冠幅和光合作用面积，使得其树木生长速度较慢生阔叶树种更快。

本研究构建了云冷杉针阔混交林中云杉、冷杉、红松、中生阔叶和慢生阔叶5个树种(组) 的枝下高和冠幅生长模型，采用了最小二乘法估计、极大似然估计和似不相关回归估计3种参数估计方法进行模型预估效果对比后，发现似不相关回归估计能够提供更加精准和稳定的参数估计值。树木胸径和树高是描述单木大小的重要变量，本研究所有树种的枝下高和冠幅均随着个体大小的增加而增加，在树木生长前期，树冠衰退速率较快，其枝下高随着单木在样地中的相对大小显著相关。当树木进入主林层后，树高和枝下高大小相对固定，树木将营养物质主要用于树冠枝条的生长，冠幅大小与林地立地质量显著正相关。尽管本研究加入了立地条件的海拔、坡度和坡向等因子，但是受样地位置相对集中的影响，仅监测到坡向对冷杉的枝下高有显著影响，在后续的研究中，可以扩大地域差异，增加立地条件差异性，以便建立适用性更广、具有更高生物学意义的枝下高和冠幅模型。