孙 圆,顾炜江,林秀云,姚睿涵,曹 林,史博文,曹福亮
(南京林业大学 a.林学院;b.南方现代林业协同创新中心,江苏 南京 210037)
落羽杉Taxodium distichum,属落叶大乔木,是生态防护、建筑用材、景观绿化等的重要树种。落羽杉自上世纪在江苏引种,目前绝大部分成林、成材,发挥了较好的生态、经济和社会效益。其中东海县李埝林场的落羽杉种源引种基地是多年来的重要研究对象[1]。前十年的研究分析显示该地引种的16 个种源的落羽杉表现出了不同的生长特性,有的种源高生长优良[2],有的种源材积生长大[3],有的种源主杆削度大[4]。又经过了近10年的经营,该实验地的落羽杉逐渐进入成过熟龄,急需在前序研究的基础上进行蓄积和材种分析,建立专用的区域树种材积表和造材表,以满足当前森林资源监测多目的多功能需要。以解析木法编制材积表的方法易受到集约经营采伐的限制,而地面激光扫描(TLS,Terrestrial Laser Scanning)能够实现立木的“实景复制”[5],通过立木点云数据能实现立木胸径、树高、上部直径和材积等测树因子的无损测量[6-8],在大量研究在表现了其时实、高效和高精度性[9-10]。特别是上部直径的提取,为材种出材率表的编制[11],为森林经营和管理提供重要了削度指标[12]。为此本研究利用地面激光扫描仪获取立木点云数据,提取研究区的落羽杉上部直径、树高、材积等测树因子,构建研究区落羽杉的一元材积模型和削度模型,编制落羽杉一元材积表和材种出材率表,为林业生产决策提供数据支撑。
研究区位于江苏省东海县李埝林场,落羽杉种源实验地中涉及16 个种源的落羽杉,共96 个小区,每个区4 行,每行6 株,株行距为2.0 m×3.3 m,采用育容器苗并定植,造林时间为1993年。前期研究于2006年得到实验区各种源生长特性结论:认为材积生长以6、8、9、10、12、13 号种源较好,各种源平均削度都较大,4、5、14、16、17 号削度性状均较好,地上生物量较大的是4、6、17、30 号种源。在此基础上本次研究着重以6 个种源(4、8、10、14、17、30 号)落羽杉纯林为研究对象。这6 个种源分别来自阿肯什州O.M.国有林、密西西比州R.F.、密西西比州5 号树、路易斯安那州N.R.、田纳西州H.C.O.,以及河南鸡公山。首先进行所有区组立木踏查,用胸径尺和测高器测定每个种源4 个小区立木测树因子,根据平均胸径选取1 株平均标准木。再采用地面激光扫描仪获取标准木及其周围邻近4 株立木的三维点云,通过点云还原立木结构。
所用扫描仪为RIEGL VZ-400i 多回波地面激光扫描仪,其扫描速率为500 000 点/s,以“井”字型布站方式进行多站扫描,如图1所示。对扫描后的数据进行了拼接、去噪、归一化、单木分割等预处理,取得每个扫描工程的5 株点云数据。最后在每个种源中取1 株平均标准木制作解析木作为验证数据。
点云数据立体还原活立木的结构现状,根据前期研究结果采用最小二乘法[9],通过散点拟合圆心,提取落羽杉立木主杆地径、胸径和上部直径,上部直径采用1 m 间隔区分段,点云的横断面切片厚度为0.1 m[12]。研究利用 FUSION/LDV 软件提取立木生长方向z 轴最大值和地面最小值的差值获得树高[13]。随后用中央断面积区分段法求得单株立木材积,其中中央断面积为点云提取的上部直径计算得到,采用1 m中央断面积区分求积法,不足1 m 区分段部分作为梢头部分,累加计算得到立木材积[14-15]。
1.3.1 立木材积表的编制与检验
一元材积表具有广泛的适用性,其以胸径为变量,在林业调查中有重要的地位。本研究选择了5 种常见的一元材积方程在DPS 17.1 软件中进行拟合,见表1。采用相关系数(R2)、平均预估精度(P)、估计值的标准误差(SEE)、平均百分标准误差(MPSE)、总相对偏差(TRB)和平均系统偏差(MSB)等评价指标,综合评价一元材积方程的拟合效果[16]。利用6 株解析木的胸径和带皮材积检验所建立的材积表适用性。
表1 一元材积表备选方程†Table 1 Alternative equations of single entry volume
1.3.2 削度方程模型与造材计算
削度方程能够描述树干形状的变化并预估不同树高处的直径,是造材计算的重要依据。研究利用立木点云数据提取的胸径(D)、树高(H)、上部直径(d),根据文献[9,12,17]的研究方法选择以下5 种削度方程(表2),在DPS 17.1 软件中建模求参。模型预测效果检验采用10 折交叉验证法,在预测模型参数已经确定的情况下将用于建模的所有单株数据随机分割成n个子样本,随机选取一个子样本作为验证数据(代入后续拟合的模型进行验证),其他n-1 个样本用来训练。验证重复10 次,每个子样本验证一次,平均n次的结果,得到模型估测值。用相关系数(R2)、残差平方和Q值优选模型作为优选指标,选出最优削度模型。
表2 削度方程模型备选方程†Table 2 Alternative equations of taper equations
造材计算根据材积规格,按照大材大用,先造大材后造小材,充分利用木材的原则进行单株木造材。将落羽杉造材数据按中华人民共和国国家标准的材种规格划分出大径材、中径材、小径材和纸浆材4 组数据[17],大、中、小径材小头直径分别不小于26.0,20.0,6.0 cm,其材长均不小于2.0 m,短小材小头直径不小于4.0 cm 且不大于12.0 cm,其材长不小于1.0 m 且不大于4.8 m。
研究区6 个种源落羽杉,每种4 个小区包括96 株立木,共有576 株踏查数据汇总如表3所示。
表3 研究区6 个种源立木数据汇总Table 3 Data collection of 6 provenances in study area
在每个种源4 个小区中对标准木进行地面激光扫描,取得每小区5 株立木点云,共120 株。从落羽杉立木点云中提取地径、胸径和树高,计算材积并按种源汇总如表4。
表4 点云测定落羽杉测树因子概况Table 4 General situation of measurement factors of bald cypress with point cloud
对编表立木胸径按5 cm 径阶制作径阶分布图,共有5~10、11~15、16~10、21~25、26~30、31~35 cm 径阶,如图2所示。全立木胸径分布呈现正态,峰值在径阶20~25 cm。对编表立木树高按2 m 树高阶制作分布图,如图3所示。全立木树高分布呈现正态,峰值在9~11 m。
图2 编表立木径阶分布Fig.2 Diameter distribution of standing trees
图3 编表立木树高分布Fig.3 Height distribution of standing trees
以实测胸径为依据检验点云提取的120 株落羽杉的胸径,如图4所示。两组数据的截距为0.141,斜率为0.99,拟合直线的R2达到0.922。表明点云提取的直径有较好的精度,可用于后续建模研究。用区分求积法得到每株立木材积,与江苏省主要树种一元材积表的数值进行对比,如图5所示。存在较为明显的偏差,总相对偏差(TRB)达到了16.08%,平均系统偏差(MSB)为33.23%,平均百分标准误差(MPSE)更是达到了41.17%以上。利用现有材积表为本实验区计算材积,整体数据偏大,不适用于实际情况。需要制定适合实验区立木状况的一元材积表。
图4 点云提取胸径与人工实测胸径散点图Fig.4 Scatter plot of DBH and measured DBH from point cloud
图5 落羽杉点云立木材积与一元材积值对比Fig.5 Comparison of the value of point cloud standing timber and the value of single entry volume
利用点云数据计算得到的材积(V)和提取的胸径(D),拟合所有备选一元材积方程。拟合结果如表5,拟合曲线如图6。从表5中可以看出,方程(2)的拟合效果最好,相关系数R²为0.969,相对偏差(TRB)在3%以内,符合《一元立木材积表编制技术规程》。因此选择方程(2)y=0.000 6x2-0.005 4x+0.034 作为建落羽杉一元材积表的方程。得到落羽杉一元材积表见附录1。
图6 一元材积模型拟合曲线Fig.6 Fit curve of single entry volume model
表5 一元材积方程拟合结果Table 5 Single entry volume models fitting result
分别为6 个种源的落羽杉立木数据拟合一元材积模型(2)。如下表6可知R2显著提高至最大0.981(4 号种源),最低0.962(8 号种源)。
表6 六个种源落羽杉一元材积y=a0+a1x+a2x2 模型拟合结果Table 6 Single entry volume model y=a0+a1x+a2x2 result for 6 provenances
利用解析木数据进行材积表检验。取得6 个种源解析木各一株,查数圆盘数据,计算伐倒木每个年龄的树干去皮材积。通过圆盘树皮数据换算成相应直径下的带皮材积。参考文献[15]方法完成材积数据。查取新建落羽杉一元材积表中的材积数据。绘制的散点图如图7。
图7 解析木数据验证一元材积表Fig.7 Analytic trees data to verify single entry volume
对模型做适用性F检验,结果如表7所示。在置信区间95%下,两组数据方差不存在显著不同(F=0.880),说明理论值与实际值差异不显著,所建立的一元材积模型符合适用性的要求。
表7 解析木材积与一元材积表材积F 检验结果Table 7 Analysis of the results of the F-test of the analytic trees volume and single entry volume volume
采用落羽杉立木点云提取的任意高度h处的带皮直径d、胸径D以及树高H值,对5 个削度方程进行拟合,结果如表8所示。其中方程(8)R2最大、残差平方和Q最小,因此选取方程(8)作为落羽杉的削度方程。
表8 各削度方程拟合结果Table 8 Fitting result of each taper equation
利用建立的削度方程推导得到累积材长方程,并计算实际造材中的大、中、小径材长度L26、L20、L6,经过积分推导得出大、中、小径材积,相关公式如表9所示。
表9 不同规格造材公式†Table 9 Volume functions with difference log rule of timbers
用以上削度方程和各材种出材材长、材积公式可进行实验区落羽杉分径阶造材的的累计材长、材积和出材率。
采用解析木数据造材,对实验区落羽杉进行造材,按径阶汇总大、中、小径材造材量,检验所建立的材种出材率表。两套数据的线性回归分析结果如表10所示,线性回归斜率均接近1,且大径材回归检验精度最好,R2都为0.99。
表10 材种出材率表检验-解析木数据vs 点云数据Table 10 Regression analysis for stem analysis versus TLS models
计算各个种源的每株材种出材量,合计每个种源样地内出材量并换算为单位面积材种出材量和材种出材率,按种源汇总如表11。在6 个种源落羽杉中,单位面积材种出材量最大的为大径材(8号种源,187.57 m3·hm-2,出材率51%),从大到小排序为8 号种源>17 号种源>14 号种源>4号种源>10 号种源>30 号种源。其次是小径材(17号种源,148.36 m3·hm-2,出材率44%),从大到小排序为17 号种源>14 号种源>8 号种源>4 号种源>10 号种源>30 号种源。
表11 各种源材种出材率核算Table 11 Accounting of the outturn rate of various provenances
已有的研究采用4 站[12]、5 站[18]或者多站[6]、单站[10]的方法完成样地立木参数估计,干型提取精度从75%到85%,本研究采用12 站井字形扫描(图1),可以保证被测立木数据完整,在不伐倒、锯断和爬树的情况下获得研究区连云港市东海县李埝林场落羽杉立木精准数据,直径提取精度达到92%(图4)并计算立木材积。本研究选择标准木对6 个种源各4 个实验区进行扫描,可以反映样地中的平均生长状态,对该方法的推广值得进一步研究,特别是在测站设置和数据提取算法上,进一步精准描述林木的三维结构。
基于点云数据分析研究区落羽杉材积与其他测树因子的关系,建立相关模型,发现研究区原有通用一元材积表不适用于该树种,其原因为:一是原有材积表非落羽杉专用;二是实验地立木经过近30 a 生长,进入成熟期需要重新评估材积生长。接着研究采用最新点云技术建立削度方程反映了立木直径随着树高变化的规律[20],反映了活立木干型状态,改变了传统测树学方法需要砍伐树木获取数据的破坏性试验[19]。沿用已有研究总结的模型形式d2=f(D,H,h)形式[12],即将模型左边统一为d2既保证了方程评价的公平性,也体现了削度方程是为了估计一定高度上直径的特性,在后期使用削度方程进行材积计算时,使用d2更为方便(表8)。
基于2006年的研究结果[4],原有16 个种源落羽杉立木生长、材积累计和干形削度均有差异,经过十多年的培育,本研究针对6 个种源完成了进一步分析。可以看出,生长表现优良的种源8号、17 号和14 号,预测出材状况较好,特别是原有削度性状良好的14 号和17 号种源表现良好(表11),这两个种源作为用材树种进行推广应用,与原有研究结果一致。原有区组实验的分析有待后续研究进行。
基于地面激光点云数据为研究区编制了落羽杉一元材积表,所建材积模型R2在0.95 以上,TRB 和MSB 均在2%以内,无较为明显的系统偏差,可用于实验区落羽杉立木材积估计。最终选定改进四参数schumacher 模型(R2=0.939)作为研究区落羽杉削度方程,该模型经过解析木数据检验,线性回归斜率均接近1,且大径材回归检验精度最好(R2=0.99)。合计每个种源单位面积材种出材量和材种出材率,单位面积材种出材量最大的为大径材(8 号种源,187.57 m3·hm-2,出材率51%),其次是小径材(17 号种源,148.36 m3·hm-2,出材率44%)。综合多性状研究分析,本结果对未来造林培育、采伐造材利用可以提供指导。