基于beta回归的迎春5号杨树树干密度混合效应模型.pdf

DOI:10.12171/j.10001522.20220450基于 beta 回归的迎春 5 号杨树树干密度混合效应模型吴新华苗铮郝元朔董利虎(东北林业大学林学院，黑龙江哈尔滨150040)摘要:【目的】探究迎春 5 号杨树在树干纵向上的木材密度影响因子和变异规律，构建迎春 5 号杨树边材、心材、树皮和树干密度混合效应 beta 回归模型，为树干生物量预测和木材材性研究提供参考。【方法】以黑龙江省尚志市 90 株迎春 5 号杨树解析木数据为基础，构建迎春 5 号杨树边材、心材、树皮和树干密度的混合效应 beta 回归模型。采用相关性分析和最优子集法筛选 beta 回归基础模型的变量；利用负二倍的对数似然值、赤池信息准则、贝叶斯信息准则、调整确定系数（Ra2）、似然比检验对收敛模型进行拟合优度的评价，利用留一交叉验证法对模型进行检验，指标为平均绝对误差（MAE）和平均绝对百分比误差；结合两种抽样方式（方案：不限定相对高；方案：限定相对高在 0.1 以下）对模型进行校正。【结果】边材、心材、树皮和树干密度不仅受到相对高的影响，还分别与胸径平均生长量、年龄、胸径密切相关，基于林木因子建立的混合效应 beta 回归模型的 Ra2分别为 0.53、0.52、0.52、0.63，MAE0.05g/cm3，与基础模型相比均提高了预测精度。边材和心材密度从树干基部往上先减小后增大，在相对高 0.2 处有拐点；树皮密度从树干基部到树梢先增大后减小，在相对高 0.6 处有拐点；树干密度沿着树干向上逐渐增大。固定相对高时，边材、心材密度都与胸径平均生长量呈负相关，树皮、树干密度分别与年龄、胸径呈负相关。在不限定相对高的情况下，沿着树干随机抽取 4 个圆盘的密度测量值来校准模型得到稳定的预测精度；限定取样高度在相对高 0.1（2.0m）以下时，对边材、心材、树皮和树干分别抽取一个圆盘（对应高度为 1.0、1.3、2.0、1.0m）的密度测量值，得到与最优抽样组合相似的预测精度。相对高、胸径平均生长量、年龄和胸径是迎春 5 号杨树木材密度的显著影响因子。【结论】beta 回归模型可对（0，1）区间的迎春 5 号杨树树干密度直接模拟，引入随机效应可提高模型的预测精度。边材、心材、树皮和树干密度在树干纵向上的变化规律不同，构建的混合效应 beta 回归模型可为迎春 5 号杨树树干生物量估算和木材性质研究奠定基础。关键词：迎春 5 号杨树；木材密度；beta 回归；广义线性混合模型中图分类号：S781.31；S792.11；S758.1文献标志码：A文章编号：10001522(2023)05006712引文格式：吴新华，苗铮，郝元朔，等.基于 beta 回归的迎春 5 号杨树树干密度混合效应模型 J.北京林业大学学报，2023，45(5)：6778.WuXinhua,MiaoZheng,HaoYuanshuo,etal.MixedeffectmodelofstemdensityofPopulus nigra P.simoniibasedonbetaregressionJ.JournalofBeijingForestryUniversity,2023,45(5):6778.Mixed effect model of stem density of Populus nigra P.simoniibased on beta regressionWuXinhuaMiaoZhengHaoYuanshuoDongLihu(SchoolofForestry,NortheastForestryUniversity,Harbin150040,Heilongjiang,China)Abstract:Objective This paper aims to explore the influencing factors and variation rules of wooddensityinthelongitudinalstemofPopulus nigraP.simonii,sobetaregressionmodelswithmixedeffectofsapwood,heartwood,barkandstemdensityofthepoplarwereconstructed,whichwasusedasareferenceforstembiomasspredictionandwoodtimberproperties.MethodMixedeffectbetaregressionmodelsforsapwood,heartwood,barkandstemdensityofP.nigra P.simonii wereestablished,whichbasedontheanalyticaldataof90treesofP.nigrax P.simoniiplantationinShangzhiCity,HeilongjiangProvinceof收稿日期:20221107修回日期:20230315基金项目:中央高校基本科研业务费专项（2572020DR03），黑龙江头雁创新团队计划项目（森林资源高效培育技术研发团队）。第一作者:吴新华。主要研究方向：木材密度研究。Email：地址：150040 黑龙江省哈尔滨市香坊区和兴路 26 号东北林业大学林学院。责任作者:董利虎，博士，教授。主要研究方向：林分生长与收获模型、生物量、碳储量。Email：地址：同上。本刊网址:http:/；http:/第45卷第5期北京林业大学学报Vol.45，No.52023年5月JOURNALOFBEIJINGFORESTRYUNIVERSITYMay，2023northeasternChina.Usingcorrelationanalysisandoptimalsubsetmethodstoscreenthevariablesofthebetaregressionbasemodel,andthegoodnessoffitoftheconvergencemodelwasevaluatedby2log-likehoodvalue,akaikeinformationcriterion,bayesianinformationcriterion,adjustedcertaintycoefficient(Ra2)andlikelihoodratiotest.Theleave-one-out-cross-validationwasusedtotestthemodel,theindexesweremeanabsoluteerror(MAE)andmeanabsoluteerrorpercentage.Twosamplingmethodswerecombined(scheme:norelativeheight;scheme:limitrelativeheightbelow0.1)tocorrectthemodel.ResultThedensitiesofsapwood,heartwood,barkandstemwerenotonlyaffectedbyrelativeheight,butalsocloselyrelatedtotheaveragegrowthofDBH,ageandDBH,respectively.Ra2ofthemixed-effectbetaregressionmodelbasedontreefactorswas0.53,0.52,0.52,0.63,respectively,andtheMAE0.05g/cm3.Sapwooddensityandheartwooddensitydecreasedfirstandthenincreasedfromthebasetothetopofthestem,withaninflectionpointatarelativeheightof0.2.Barkdensityfirstincreasedandthendecreasedfromthebaseofthestemtothetopofthetree,andtherewasaninflectionpointattherelativeheightof0.6.Thestemdensityincreasedgraduallyalongthestem.Whenfixedrelativeheight,thedensitiesofsapwoodandheartwoodwerebothnegativelycorrelatedwiththeaveragegrowthofDBH.ThedensitiesofbarkandstemwerenegativelycorrelatedwithageandDBH,respectively.Withoutlimitingtherelativeheight,thewooddensityvaluecorrespondingtotheheightof4discsrandomlysampledalongthestemwascalibratedtoobtainstablepredictionaccuracy.Whenthesamplingheightwaslimitedto0.1(2.0m)orless,therewaslittledifferenceinthepredictionaccuracybetweentheoptimalsamplingcombinationandthedensityvalues(1.0,1.3,2.0,1.0m,respectively)ofsapwood,heartwood,barkandstematadischeight.Relativeheight,averagegrowthof DBH,age and DBH were significant influencing factors of wood density of P.nigra P.simonii.ConclusionThebetaregressionmodelcandirectlysimulatethestemdensityofP.nigra P.simoniiinthe(0,1)interval,andtherandomeffectcanimprovethepredictionaccuracyofthemodel.Thelongitudinalvariationsofsapwood,heartwood,barkandstemdensityaredifferent.Theconstructedmixed-effectbetaregressionmodelcanlayafoundationforbiomassestimationandwoodpropertystudyofP.nigra P.simonii.Key words:Populus nigra P.simonii;wooddensity;betaregression;generalizedlinearmixedmodel木材密度又称为木材基本密度，通常用绝干材质量除以生材体积来表示。木材密度与碳储量、木材机械强度、水分运输、耐旱性等密切相关12，在许多研究利用树干木材密度乘以树干体积来计算树干生物量34，进而估算森林碳储量5。碳储量与生态环境变化息息相关，所以木材密度对研究森林生态系统固碳功能具有重要作用6。木材密度也是林木生长率和死亡率的重要预测因子，可以反映林木的生长和性状表现7，通常低木材密度代表高生长量810。准确估算木材密度不仅有助于阐明地上生物量的空间分配格局，还能反映树干内部的机械稳定性和生长速率的变化1112。作为森林碳汇计量中的重要生物量估算因子13，木材密度主要受树种、相对高、树龄、生长量、胸径和环境条件等因素的影响1417。有研究表明杨树（Populusspp.）的木材密度沿树干自下往上逐渐增大18，而落叶松（Larix gmelinii）却呈现出相反的规律。樟子松（Pinus sylvestris）、日本花柏（Chamaecyparispisifera）等树种木材密度沿树干自下往上先减小后增大1920，这表明不同树种间的木材密度变化规律具有较大差异。Iida 等21的研究表明树高一定时，较小的胸径具有较高的木材密度。木材密度还与年轮宽度，尤其与晚材比例密切相关，在环孔材中更加明显22。研究表明干旱条件会抑制树木生长进而增加晚材比例，获得较大的木材密度，进一步说明低含水率导致较高的木材密度23。目前对松属（Pinus）、杉科（Taxodiaceae）以及硬阔树种木材密度的研究较多19,2425，暂未见迎春 5 号杨树（P.nigra P.simonii）木材密度的相关报道。各研究显示不同树种之间木材密度有所差异，同一树种内部木材密度也不尽相同，例如边材、心材和树皮木材密度；径向和纵向木材密度；幼树和成熟树木材密度26。准确把握树干各部分的木材密度变化规律，合理预测木材密度，有利于生物量的精确预估。早期的研究者常利用线性和非线性回归来构建木材密度模型，种类有多项式、双曲线、指数方程等27。但木材密度数据位于（0，1）区间内，在数据结构上属于连续有界数据，利用上述模型对该类数据68北京林业大学学报第45卷进行拟合时，容易产生超过其上下限的值，降低预测精度28。此外，（0，1）区间的数据存在分布不对称的问题，无法基于正态性假设来进行推断预测。广义beta 回归模型非常灵活，能够直接应用于连续有界数据，并且可以适应“J”形、“L”形、“单峰”、“U”形和矩形的各种类型的分布数据29，解决了数据不对称分布问题。另外，木材密度在不同树木、同一树木的不同位置间具有较大的变异性，且数据一般来源于固定样地或样木，数据之间存在一定的相关性，并不满足传统回归模型中独立同分布的假设条件，故利用混合效应模型对木材密度进行建模，通过随机效应方差协方差结构来表现个体间的差异和个体内的自相关性，提高模型预测精度3032。迎春 5 号杨树是小叶杨（P.simonii）和黑杨（P.nigra）杂交而成的优良无性系。该品种自 1982 年从林科院引入牡丹江林区，是东北地区主要造林树种之一，同时也是培育速生丰产林及短轮伐期工业原料林的优良品种。基于经济效益，前人对于迎春5 号杨树的研究主要集中在造林方法、种植密度和经营技术方面，暂未见其木材密度以及生物量的研究。本研究从生态效益出发，构建迎春 5 号杨树木材密度模型，进而估算其树干生物量，为杨树人工林碳储量和木材材性研究奠定基础，具有生态学意义。1研究区概况与研究方法 1.1 研究区概况研究区位于黑龙江省尚志市的尚志林场（450008452328N，12739321293509E）、元宝林场（450826451826N，12810171283307E）、老街基林场（440906450401N，12744091280506E）和苇河林场（440109445201N，12810521284704E）。各林场隶属于尚志国有林场管理局，地处长白山脉张广才岭西坡，地形以山地丘陵为主，地势西低东高，土壤以暗棕壤为主。环境气候条件相似，均为中温带大陆性季风气候，冬季寒冷干燥，夏季温暖湿润，降水量充沛。森林资源丰富，尚志国有林场森林总面积为152737hm2，其中人工林面积 49168hm2，占比32.19%。植被属长白山植物区系，人工林主要树种为长白落叶松（Larix olgensis）、樟子松、红松（Pinuskoraiensis）和杨树。2020 年 9 月在尚志市上述林场选择不同年龄、不同林分密度、不同海拔的迎春 5 号杨树人工林设置标准地共 18 块，其中尚志林场 11 块，元宝林场3 块，老街基林场 1 块，苇河林场 3 块。标准地的形状为矩形。由于林分疏密程度不同，所设标准地面积 0.060.10hm2不等。各标准地所在林场的相关立地因子如表 1 所示。在每块标准地中进行每木检尺，按照等断面积标准木法将每块标准地内的林木划分为 5 级，径阶从大到小依次为优势木、一级木、二级木、三级木和劣势木。以每级木的平均胸径为标准，在每块标准地附近选取 5 棵胸径大小误差为5%左右的林木作为解析木，对其进行树干解析，共获得 90 棵解析样木。表1各林场相关的立地因子Tab.1Sitefactorsrelatedtoeachforestfarm林场Forestfarm地形Terrain海拔Altitude/m坡度Slope坡向Aspect坡位Slopeposition尚志林场ShangzhiForestFarm山坡Hillside20131905、515阳、阴Sunny,shady中下Lower-middle元宝林场YuanbaoForestFarm平地Flat20624505平Flat老街基林场LaojiejiForestFarm平地Flat30005平Flat苇河林场WeiheForestFarm山坡Hillside26032405、510阳、阴Sunny,shady中下Lower-middle在进行树干解析前，要在树干上标记北向和1.3m 处，并测量胸径。解析木伐倒后，测量其树高和年龄等因子。从解析样木根颈（0m）处每隔 1.0m截取 2cm 厚的圆盘样品，截取过程中避免节子，胸径（1.3m）处单独取盘，直至距离树梢不足 1.0m 为止。按高度对圆盘进行编号，在非工作面标注编号和北向，将不同解析样木的圆盘分别装袋带回实验室。各圆盘以髓心为定点，取适当大小的楔形或全部圆盘作为该圆盘的样品。由于迎春 5 号杨树的边材、心材和树皮颜色差异明显，因此将样品分为边材、心材和树皮 3 部分。若样品没有心材，则将其分为边材和树皮两部分。利用排水法分别测量边材、心材和树皮样品入水前后的质量，利用水的密度分别计算排出水的体积，即为边材、心材、树皮样品的生材体积。最后将样品置于 85 烘箱 24h，测量样品的绝干质量，木材密度表示为每单位生材体积的干质量，用绝干材质量除以生材体积来计算33。迎春 5 号杨树人工林主要林分因子和解析木因子信息见表 2，各计算公式如下。第5期吴新华等：基于 beta 回归的迎春 5 号杨树树干密度混合效应模型69V=(m2m1)/0(1)m=m/V=m0/(m2m1)(2)m1、m2m式中：分别代表样品入水前后的质量（g）；0代表水的密度（g/cm3）；代表木材密度（g/cm3）；m 代表样品绝干材质量（g）；V 代表样品的生材体积（cm3）。1.2 研究方法1.2.1beta 回归基础模型beta 回归模型可对分布在（0，1）区间上的连续变量进行模拟，模型较为灵活，故本研究将其作为预测迎春 5 号杨树纵向上边材、心材、树皮和树干密度的基础模型。模型的因变量为边材密度、心材密度、树皮密度和树干密度。模型预选自变量为相对高、树高、胸径、树龄、胸径平均生长量、林分平均胸径、林分平均高。具体建模过程在 SAS 软件 glimmx 模块中进行。beta 分布的概率密度函数34如下。f(y;,)=()()(1)y1(1y)(1)1(3)0 y 1()E(y)=，0 1 var(y)=(1)1+式中：，为伽马函数；参数 表示 y 的期望值，；为分散度参数且为正数，越大，y 的方差越小。可以表示为y 的方差。Yiig(i)=YiXii=1，2，nXiYi(0，1)Yi betai，(1i)g(i)=XTi YiXi与其他的广义线性模型相似，beta 回归模型的响应变量的期望值是由一个链接函数构造的，是关于响应变量的线性组合。假设有协变量，且 n 为正整数，的响应变量且，则链接函数（为回归系数向量10）。链接函数的引入是为了将响应变量与协变量之间的非线性关YiiXig(i)=lni1i=XTi系转换为 的期望与间的线性关系，在不改变数据原始结构的基础上直接获得所需的拟合值。常见的链接函数有 logit、log 和 probit 等28。因为 logit 链接函数中的回归系数可以转化为优势比，能够更好地解释响应变量与协变量之间的关系10,35，故本研究中运用 logit 链接函数，即。采用极大似然法进行参数估计。1.2.2混合效应 beta 回归模型在 beta 回归模型的基础上加入随机效应，构成混合效应 beta 回归模型。由于数据中样地间差异较小，因此本研究未考虑样地水平的随机效应，仅考虑样木水平的随机效应。单木水平的混合效应 beta 回归模型包括固定效应参数和随机效应参数，其中固定效应参数是对总体均值变化的反映，随机效应参数则是对个体变异程度的反映，在本研究中具体表现为不同样木间的木材密度差异。在线性混合模型中，通常需要假设随机效应和方差服从正态分布34，但本研究中的木材密度测量值位于（0，1）区间且有较大的变异性，并不满足上述条件。本研究采用混合效应 beta 回归模型来对迎春 5 号杨树边材、心材、树皮和树干密度进行建模，模型形式如下。g(i)=lni1i=XTi+ZTibi(4)bi N(0，D)DXiZibibi式中：，表示随机效应的方差协方差矩阵。为固定效应的设计矩阵；为固定参数向量；为随机效应的设计矩阵；为随机参数向量；模型通过对积分获取最大似然值，进而得到参数估计值36。构建混合效应 beta 回归模型时，首先在不同位置加入随机效应参数比较拟合优度；其次确定随机表2林分因子和解析木因子基本信息统计表Tab.2Statisticaltableforbasicinformationofstandfactorsandanalyticalwoodfactors因子Factor最小值Min.value最大值Max.value平均值MeanSD林分因子Standfactor林分年龄/aStandage/year1026185平均胸径AverageDBH/cm13.330.620.65.3平均树高Averagetreeheight/m11.926.220.34.1林分密度/（株hm2）Standdensity/(treeha1)183933572267树木因子Treefactor树木年龄/aTreeage/year1026185胸径DBH/cm6.237.521.66.4树高Treeheight/m8.928.921.14.2边材密度Sapwooddensity/(gcm3)0.270.530.360.03心材密度Heartwooddensity/(gcm3)0.260.500.330.03树皮密度Barkdensity/(gcm3)0.200.550.370.07树干密度Stemdensity/(gcm3)0.270.510.360.0370北京林业大学学报第45卷Ri效应方差协方差矩阵（D），林业上常用的 3 种方差协方差结构分别为：复合对称（compoundsymmetry，CS）、对角矩阵（diagonalmatrix，UN（1）、正定矩阵（generalpositive-definitematrix，UN）；最后确定组内方差协方差结构（）。Ri=2G0.5iiG0.5i(5)2Gii式中：为模型的残差方差，为描述方差异质性的对角矩阵，为树木内自相关矩阵。1.2.3模型评价与检验用全部数据拟合模型，模型拟合效果的评价指标为负二倍的对数似然值（2loglikelihoodvalue，2lnL）、赤池信息准则（Akaikeinformationcriterion，AIC）、贝叶斯信息准则（Bayesianinformationcriterion，BIC）和调整确定系数（adjustedcertaintycoefficient，Ra2），前三者数值越小，后者数值越大，模型拟合效果越好。利用留一交叉验证法（leave-one-out-cross-validation，LOOCV）对模型进行检验，即每次抽取1 条检验数据，其他数据作为建模数据重新拟合，一直重复直到所有数据均被验证。检验指标为平均绝对误差（meanabsoluteerror，MAE）和平均绝对百分比误差（meanabsolutepercentageerror，MAPE），二者数值越小模型越好。模型参数个数不同，拟合效果也不同，利用似然比检验（likelihoodratiotest，LRT）比较不同参数个数对模型拟合效果的影响，若LRT 的显著性 P0.05，表示两模型间有显著性差异。各指标的公式如下。AIC=2ln L+2(6)BIC=2ln L+ln n(7)R2a=1n1nni=1(YibYi)2ni=1(YiY)2(8)ma=ni=1?YibYi,in?(9)P=1nni=1?YibYi,i?Yi100%(10)RL=2ln(L2L1)=2(ln L2ln L1)(11)maPYibYiYbYi,iL1L2式中：AIC 为赤池信息准则；BIC 为贝叶斯信息准则；L 为似然函数；2lnL 为负二倍的对数似然值；为平均绝对误差；为平均绝对百分比误差；RL为似然比；为模型参数个数；n 为样本量；为密度实际值；为密度预测值；为密度测量值的平均值；为刀切法之后的密度预测值；、分别为一般模型和限制模型的似然值。1.2.4模型预测在混合效应模型中，需要对数据进行二次抽样来计算随机效应参数估计值对模型进行校正，以获取更高的预测精度。本研究采用林业上常用的最佳线性无偏估计3738来计算随机效应参数。bbk=bDbZTi(bZkbDbZTk+bR)1b ek(12)bbkbDbZkbRb ek式中：为随机参数估计值，为随机效应参数的方差协方差矩阵，为设计矩阵，为组内方差协方差矩阵，为实际值减去固定效应参数计算的预测值。考虑到实际操作中的成本问题，本研究设计了两种抽样方案对模型进行预测。方案：不限定相对高时，随机抽取任意位置的圆盘木材密度值，获取预测精度稳定时的抽样数量。方案：限定相对高 0.1（2.0m）以下时，从每棵树中抽取 14 个圆盘的木材密度值，获取不同抽样数量时模型预测精度最高的抽样组合并比较。2结果与分析 2.1 beta 回归基础模型对圆盘数据进行分析后发现：与树高 H 相比，相对高 Hr更能反映木材密度纵向上的变化规律，因此将 Hr或者 Hr2作为模型固定变量。根据其他预选变量间的相关系数（图 1 中显示的相关系数值代表变量间显著相关，不显示数值表明变量间无显著相关性）以及最优子集筛选后发现，胸径平均生长量 Dt、年龄 t、胸径 DBH与木材密度关系密切。通过比较2lnL、AIC、BIC、Ra2，并对参数不同的模型进行似然比检验，最终确定预测迎春 5 号杨树的边材密度基础模型的变量依次为 Hr2、Hr、Dt；预测心材密度基础模型的变量依次为 Hr2、Dt、Hr；预测树皮密度基础模型的变量依次为 Hr、Hr2、t；预测树干密度基础模型的变量为 Hr2、DBH、Hr。基础模型通式如下。E(Y)=1/1+exp(0+1x1+2x2+3x3)(13)x1、x2、x30代表截距，1、23式中：Y 分别代表边材密度、心材密度、树皮密度和树干密度；分别代表各个密度模型中对应的变量；、分别代表自变量前的系数。2.2 混合效应 beta 回归模型的建立与最优模型的选取为了反映不同树木间的差异，在不同位置加入随机效应参数，比较模型的拟合优度，进而选择最优的单木水平混合效应 beta 回归模型来探究不同树木间的边材、心材、树皮和树干密度的变化规律。模型拟合结果如表 3 所示。第5期吴新华等：基于 beta 回归的迎春 5 号杨树树干密度混合效应模型71由表 3 可知：边材、心材、树皮和树干密度模型分别有 2、4、4、5 种收敛情况。随着随机效应的加入，混合模型2lnL、AIC、BIC 的值较基础模型均有所减小，Ra2增大，说明混合效应模型提高了预测精度。另外，随机参数的个数对模型拟合优度有所影响，除了边材密度混合效应模型随机参数个数相同之外，心材密度、树皮密度和树干密度的混合效应模型均显示出两个随机效应参数的拟合效果比一个好。根据表 3，对于边材密度，只有一个随机效应参数时模型才收敛，对比模型评价指标，将模型 2 作为边材密度混合效应 beta 回归最优模型。对于心材、树皮和树干密度，当有一个随机参数时，各密度模型中皆为模型 2 效果最好，当有两个随机参数时，各密度模型中模型 4 效果最好。对各密度模型中的模型2 和模型 4 进行似然比检验，结果表明两模型间均有显著性差异。因此将各密度的模型 4 分别作为心材、树皮和树干密度混合效应 beta 回归最优模型，其固定参数估计值、随机效应方差协方差结构如表 4所示。UN(1)21212122比较 UN、CS 和的效果（以心材密度混合效应模型为例），3 个矩阵依次代入模型得到的2lnL分别为6383.05、6360.02、6381.04，显示 UN的效果最好，故将其作为随机效应的方差协方差矩阵，表示为（以两个随机效应参数为例）。又测试了林业上常用的两种自相关结构：AR（1）和ARMA（1，1），自相关结构的加入使得模型无法收敛，故本研究中没有进一步考虑纵向各圆盘密度间的自相关性。123123由表 4 可看出：边材、心材、树皮密度模型中的系数、和的符号均为正、负、负，树干密度模型中的系数、和的符号为正、负、正。对应各r1.00.500.51.0sDtDtHDHDHDgDgDBHDBHttr1.00.500.51.0bDtDtHDHDHDgDgDBHDBHttr1.00.500.51.0hDtDtHDHDHDgDgDBHDBHttr1.00.500.51.0wDtDtHDHDHDgDgDBHDBHtt0.060.060.110.100.090.290.350.160.180.370.300.130.080.110.110.310.390.170.070.320.390.180.070.080.240.770.630.820.860.110.420.230.750.610.810.860.640.720.690.900.690.900.120.090.660.730.700.900.690.910.420.070.430.260.760.630.820.910.860.660.730.710.910.700.430.250.750.610.810.850.640.700.690.900.690.90s为边材密度，h为心材密度，b为树皮密度，w为树干密度，t 为树龄，H 为树高，Dt为胸径平均生长量，Dg为林分平均胸径，DBH为胸径，HD为林分平均高。蓝色代表负相关，红色代表正相关，P0.05。sissapwooddensity,hisheartwooddensity,bisbarkdensity,wisstemdensity,tistreeage,Histreeheight,DtistheaveragegrowthofDBH,DgistheaverageDBHofstand,DBHisthediameteratbreastheight,andHDistheaverageheightofstand.Bluerepresentsnegativecorrelation,andredrepresentspositivecorrelation,P0.05.图1变量间的相关关系Fig.1Correlationsamongvariables72北京林业大学学报第45卷密度模型中的自变量顺序（2.1 中有所叙述）可知 Dt、t 和 DBH前的系数都为负，说明随着 Dt、t 和 DBH的增大，边材、心材、树皮和树干密度分别都减小。边材、心材、树皮密度模型中的 Hr和 Hr2前的系数一正一负，说明对应部位的木材密度在树干上存在最小值。树干密度模型中 Hr和 Hr2前的系数均为正，说明树干密度随着相对高的增加而增加。以建模数据中的 Dt、t、DBH数据的最大值、平均值和最小值作为不同梯度，根据表 4 中的各木材密度模型的固定效应参数，绘制迎春 5 号杨树边材、心材、树皮和树干密度与各变量间的关系曲线图（图 2）。图 2 中横坐标表示相对高，数字 0 和 1.0 分别代表树干基部和树梢。图 2 显示：迎春 5 号杨树边材密度从树干基部到树梢先缓慢减小后增大，在相对表3边材、心材、树皮和树干密度混合效应 beta 回归模型的拟合优度比较Tab.3Goodnessoffitcomparisonofbetaregressionmodelsformixedeffectsofsapwood,heartwood,barkandstemdensity部位Part模型编号ModelNo.随机效应Randomeffect评价指标EvaluationindexHrHr2DttDBH2lnLAICBICRa2RLP边材Sapwood07988.277978.277950.600.3118316.718304.718289.710.5028417.798405.798390.790.53心材Heartwood06012.416002.415976.120.1916237.856225.856210.850.4126324.636312.636297.630.4636291.246279.246264.240.4446383.056367.056347.050.5258.420.001树皮Bark04887.394877.394850.110.2915069.365057.365042.360.4425106.945094.945079.940.4635030.855018.855003.850.4245163.775147.775127.770.5256.830.001树干Stem07377.797367.797340.530.3617566.357554.357539.350.5027736.637724.637709.630.5637651.707639.707624.700.5347840.397824.397804.390.6357822.197806.197786.190.62103.760.001注：代表在该变量上添加随机效应，代表模型固定效应的自变量。Hr代表相对高，Hr2代表相对高的平方，Dt代表胸径平均生长量，t代表树龄，DBH代表胸径，2lnL代表负二倍的对数似然值，AIC代表赤池信息准则，BIC代表贝叶斯信息准则，Ra2代表调整确定系数，RL代表似然比。Notes:standsforaddingrandomeffectstothisvariable,standsfortheindependentvariablesofthemodelfixedeffect,Hrstandsfortherelativeheight,Hr2standsforthesquareofrelativeheight,DtstandsfortheaveragegrowthofDBH,tstandsforthetreeage,DBHstandsfordiameteratbreastheight,2lnLstandsforthe2loglikelihoodvalue,AICstandsfortheakaikeinformationcriterion,BICstandsfortheBayesianinformationcriterion,Ra2standsfortheadjusteddeterminationcoefficient,RLstandsforthelikelihoodratio.表4边材、心材、树皮和树干密度最优混合效应模型的固定效应参数估计值、随机效应方差协方差结构Tab.4Estimatedvaluesoffixedeffectparameters,randomeffectvariancecovariancestructureforoptimalmixed-effectsmodelofsapwood,heartwood,barkandstemdensity参数Parameter边材Sapwoo