基于机载激光雷达冠层高度模型的小班区划

熊 昊，庞 勇＊，荚 文，李春干

(1.中国林业科学研究院资源信息研究所，国家林业和草原局林业遥感与信息技术重点实验室，北京 100091；2.广西大学林学院，广西 南宁 530004)

小班是内部特征基本一致，与相邻地段有明显区别而需要采取相同经营措施的森林地块或小区，是森林资源统计、经营和管理的基本单位[1]。小班区划通常是依据高分辨率遥感影像目视判读或实地勾绘，工作量巨大，且会因调查员的知识水平、经验、调查线路、观察视角等主观因素差异产生结果不一致的问题[2-3]。因此，对小班自动区划进行研究十分必要。在小班的划分条件[1]中，权属、森林类别等因子是政策规定或造林时已经确定的，这些因子通常不因林分状态的变化而改变，且可以从上一期林相图中提取。而优势树种组成、龄级等因子可以通过遥感方法提取。

国内已有一些研究采用遥感影像来进行自动或半自动的小班区划。这些研究通常采用Landsat、SPOT5、QuickBird、ALOS 等高分辨率卫星遥感影像作为辅助数据，使用多尺度分割和分类的方法获得小班界线[3-6]。但卫星影像受到空间分辨率的约束，对小班区划仍具有一定的局限性。国外有一些学者采用航空遥感影像进行小班区划[7-8]，但多限于区分内部较均质的林分，没有对我国森林调查规划规程中规定的优势树种、小班面积等因子的严格界定。并且影像数据中不包含林分结构信息，可能会有树种相同，但树高、密度等不同的林分被划分为同一小班，同时影像也可能有云和阴影等干扰因素。

近年来机载激光雷达（ALS）技术在林业中的应用日益广泛[9]。ALS 数据是由搭载在有人机或无人机平台上的激光雷达传感器扫描获取的，适用于森林结构观测[10]。由ALS 数据生成的森林冠层高度模型可以为小班区划提供准确的位置、林分平均高、林分密度、郁闭度等信息[11-12]。通过对CHM的分割可以将森林参数不同的地块区分开，得到符合林区实际情况的小班边界。Mustonen 等[13]比较了采用CHM 和RGB 航空遥感影像和综合两种数据的林分区划，发现基于CHM 的区划的林分比基于RGB 航空遥感影像和基于两种数据区划的林分更均匀。Koch 等[11]利用20 m 空间分辨率的低密度ALS 数据进行小班区划，采用特征提取、创建和基于栅格的分类方法。Wu 等[14]从点云中提取树的大小、林分密度和树种特征，使用非监督均值漂移算法和区域生长法进行分割。Dechesne 等[15]从点云和多光谱影像中计算并筛选出了94 个特征，并使用随机森林方法进行树种监督分类得到林分区划结果。Pukkala 等[16-18]利用16 m 的ALS 数据和森林资源调查数据，使用元胞自动机方法，通过调整栅格内部属性参数来进行林分区划。而国内现有研究中，仅有少数研究利用ALS 数据在ArboLiDAR软件中进行了林分高度分割[19]，并用于反演云杉次生纯林的林分平均高。Jia 等[20]讨论了不同形状和面积参数对元胞自动机林分区划的影响。

已有的研究中体现了ALS 数据用于林分区划的潜力，但这些研究没有针对小班区划需求，并且未充分考虑树种信息。机载高光谱影像具有较高的空间和光谱分辨率，大量研究表明[21-25]，机载高光谱数据可以产生高精度的森林树种分类结果，为小班区划提供准确的树种组成信息。因此，本研究利用机载激光雷达数据和高光谱影像树种分类结果进行小班区划，并与人工区划小班进行对比，综合探究在高光谱影像树种信息的辅助下，ALS 数据生成的高精度CHM 在小班区划中的应用潜力。

1 研究区概况及数据

1.1 研究区概况

研究区位于黑龙江省佳木斯市桦南县东北部的孟家岗林场，地理坐标为130°32′～130°52′ E，46°20′～46°30′ N。地处完达山西麓余脉，以低山丘陵为主，坡度较为平缓，坡度在10°～20°之间。地势东北高，西南低，平均海拔约为250 m。孟家岗林场人工林占林地面积的76.7%，其中落叶松（Larix olgensisHenry）、樟子松（Pinus sylvestrisvarmongolicaLitvin.）、红松（Pinus koraiensisSieb.et Zucc.）约占80%。天然林主要分布在林场的东北方向，占全部森林面积27.5%，主要树种包括白桦（Betula platyphyllaSuk.）、柞树（Quercus mongolicaFisch.Ex Ledeb）和椴树（Tilia tuanSzyszyl）。

1.2 数据

1.2.1 遥感数据 采用中国林业科学研究院的机载综合遥感观测平台CAF-LiCHy[26]于2017 年5 月31 日至6 月15 日采集的LiDAR 数据，航摄平台为运-5 多用途飞机，绝对飞行高度为1 000 m。CAF-LiCHy 机载系统包含全波形机载激光雷达（LMS-Q680i）、机载推扫式高光谱扫描仪（AISA EagleⅡ）和高分辨率CCD 相机，3 个传感器共用一套POS 位置和姿态系统。飞行时航向和旁向交并比分别为70%和50%，获取的数字正射影像空间分辨率为0.1 m。LMS-Q680i 激光传感器是奥地利Riegl 公司的全波形激光雷达系统，工作波长为1 550 nm，激光束发散角为0.5 mrad，波形数据的记录间隔为1 ns，最大脉冲重复频率为400 kHz。激光雷达点云数据的平均密度为12 pts·m-2，航带间点云误差为水平0.1 m，垂直0.12 m。高光谱AISA Eagle II 传感器光谱范围400～970 nm，光谱分辨率为3.3 nm，空间分辨率为0.5 m。CAFLiCHy 机载系统同步获取的激光雷达数据、数字正射影像和高光谱影像的空间位置的一致性精度优于1 m[27]。

1.2.2 树种分类图 树种类别图由CAF-LiCHy 机载系统同步获取的高光谱影像监督分类获得。对高光谱影像的每个波段使用灰度共生矩阵计算了均值、均匀性、异质性、相关性、方差、二阶矩、熵、对比度8 种纹理特征，对纹理特征用最小噪声分离法降维。使用提取的纹理信息和光谱信息进行SVM 监督分类，分为落叶松、红松、樟子松、云杉、阔叶树、耕地和其他7 类[24-25]。经2017、2019 和2020 年采集的679 个样地调查数据验证，总体分类精度为91.28%，Kappa 系数为0.88。

1.2.3 样地实测数据 2017 年对孟家岗林场70 块固定样地进行每木检尺，测量了样地中每株树的胸径、树高和冠幅，计算每个样地的算术平均胸径和算术平均高。采用差分GNSS 对样地进行精确定位，定位精度优于1 m。70 个样地中包括落叶松样地40 个、樟子松样地11 个和红松样地19 个。

1.2.4 人工区划小班数据 2016 年以卫星影像为底图，通过目视解译得到林场林相图，与航空遥感数据存在时间和空间位置上的差异。因此，以2017年机载高光谱影像为底图，对2016 年林相图进行修正，得到2017 年人工区划小班数据。修正调整后的小班界线更符合实际情况，且保证了与机载高光谱和激光雷达数据空间位置的一致性。冠层高度模型、树种类别图、小班边界和样地分布见图1。

图1 孟家岗林场冠层高度模型、树种类别图及小班区划图Fig.1 Canopy height model,tree species map and related data for sub-compartments delineation in Mengjiagang forest farm

采用孟家岗林场2017—2019 年共有149 个伐区调查数据，面积为1～17 hm2。由于林分生长慢，对小班区划而言，2017—2019 年的林分变化可以忽略不计，故采用3 个年度的主伐作业小班作为对照。

在森林资源调查监测工作中，参考对象主要有人工区划的小班图和比用于自动区划的遥感数据分辨率更高的数据勾绘得到的小班图，其中勾绘得到的小班图误差更少[4]。本研究基于0.1 m 空间分辨率DOM 屏幕勾绘了100 个小班作为参考小班。

2 研究方法

使用机载高光谱影像分类得到的树种图作为辅助数据，首先对1 m 空间分辨率的CHM 数据进行多尺度分割，得到过分割结果图；然后对1 m 空间分辨率的CHM 数据均值下采样到5 m，并使用对称近邻均值滤波（SNN）平滑，再对其进行多尺度分割，得到欠分割结果图；最后将过分割和欠分割结果进行融合优化，得到最终小班区划结果。总流程图如图2 所示。

图2 基于机载激光雷达冠层高度模型和树种类别的小班自动区划流程Fig.2 The flowchart of automatic delineation of forest stands based on CHM from airborne LiDAR data and tree species

2.1 LiDAR 数据预处理

采用Lastools 软件对LiDAR 点云数据进行预处理。去除噪声点，进行点云分类，将分类后的地面点与非地面点分别采用TIN 内插算法与最大高程插值方法生成空间分辨率为1 m 的数字高程模型（DEM）和数字表面模型（DSM），用DSM与DEM 相减获得CHM 数据，其空间分辨率为1 m[28]。

2.2 CHM 平滑滤波

CHM 空间分辨率小于研究区树冠的平均冠幅，在图像分割中将产生大量过分割的斑块，为此，对1 m 空间分辨率的CHM 数据进行适当的重采样和平滑滤波，对同一株树的多个像元取平均值，以更好地显示这株树的信息。在林分尺度上，对CHM 数据适当的平滑不会对小班边界产生大的影响。为减少像元重采样和平滑对小班边界的影响，确保边界的突显性，采用具有边缘保持特性的平滑滤波器。

对KNN 均值、最小方差均值、Sigma 平滑和SNN4 种滤波器的5 m 空间分辨率CHM 数据平滑效果进行比较分析，综合考虑异质区域间边界的保持效果和同质区域内的平滑效果，选用SNN 滤波器进行5 m 分辨率CHM 的平滑处理。

SNN 滤波的核心思想是：比较窗口范围内对称像元的值，选取每对像元中与目标像元值更相近的值，对窗口内选取的所有像元取均值作为目标像元值，并保留与目标像元差异大的像元值不变。这样在保持边缘的同时，又可以以较低的计算量对图像进行平滑[29]。根据地面调查数据，研究区成熟林分的平均冠幅约为5 m。比较3 × 3、5 × 5、7 ×7 和9 × 9 像素窗口的滤波效果，发现当窗口过小时会使平滑效果变差，最终区划结果中产生过多面积较小的小班。窗口过大时不同小班的边界更加模糊，难以区分不同小班，最终区划结果中有很多差异较大的应区划为不同小班的区域未分割开，窗口大小为5 × 5 时滤波结果最优。

2.3 多尺度分割

多尺度分割算法是从单像元对象开始的区域合并过程，合并规则基于局部光谱均匀性，即相邻图像对象之间的相似性。光谱均匀性标准定义为光谱值的标准偏差之和，每次合并时方差增加最小的对象被合并，当最小增量高于用户定义的阈值(尺度参数)时，分割过程结束[30]。本研究中的过分割即分割出大量面积较小的对象，对象小于小班尺度。欠分割即分割出数量较少、面积较大的对象，对象接近小班尺度。其中面积很小的与周围区域差异较大的区域没有被分割出来，因此称为欠分割。

将人工区划小班叠加到树种分类图上，统计每个小班各树种像元数占小班总像元数的比例，即树种组成。选取树种成数最大的优势树种作为该小班的树种。融合相同树种的区域形成树种矢量图。在eCognition 软件中用CHM 数据作为底图，加入树种矢量图作为辅助数据，在同类树种区域内对CHM 进行多尺度分割得到林分高度和密度一致的小班斑块。分割方案中尺度参数设定为10～300，形状参数和紧致度参数设定为0.1～0.9，进行了多组试验，最终确定了分割效果最优的参数。

2.4 分割结果优化

对CHM 数据进行多尺度分割的过程中，采用了多组参数进行试验，结果显示为1 m 和5 m 空间分辨率数据的分割结果最具互补性。5 m 空间分辨率CHM 的欠分割结果显示了主要的林分区域，但其轮廓和林分边缘的树冠轮廓并不完全相符。1 m分辨率CHM 的过分割结果虽然形成了过多的斑块，但其轮廓和林分边缘的树冠轮廓更相符。因此需要使用过分割结果来优化欠分割结果。

本研究使用ArcGIS 软件中的update 工具实现多尺度分割结果的融合。其原理是对输入图层和用于更新的图层进行几何相交的计算，其中输入图层被更新图层覆盖部分的属性用更新图层的属性代替。在更新时可以调整容差值，即输出图层边界的X,Y坐标相对于输入图层边界的X,Y坐标可移动的距离。将欠分割结果作为更新图层，通过调整两种尺度的分割结果间的容差值，经过多次试验，确定容差值为3 m 时优化效果最佳。

2.5 区划结果精度验证

区划结果采用形状因子和定量因子2 种方式进行精度验证。首先评估林分边界的准确程度，采用人工区划小班、主伐作业小班和0.1 m 空间分辨率DOM 屏幕勾绘的小班作为参考。根据自动区划小班与参考小班的空间位置相对应的原则，采用最终测量精度准则，逐一测量自动区划小班与参考小班的特征。最终测量精度准则评价的特征包括圆度（RO），紧致度（CO），形状指数（SI），最小包络圆短半径（RE），椭圆度（EF）和形状因子（P2A）[4]。

同时，使用交并比指标[31-32]，比较自动区划小班和参考小班的重合情况。计算人工区划小班与其对应的自动区划小班的交集的面积与并集面积的比值[31]，即：

其次，采用样地平均胸径、平均树高和5 m ×5 m CHM 的冠层平均高数据计算可解释性方差精度，评价小班内部的均质性及不同小班之间的异质性[16-17,20]。即计算小班内部观测值或冠层平均高的方差之和SSwithin，各样地值与全研究区均值的方差之和SStotal的可解释性方差R2，公式如下：

使用样地数据计算时，k为小班数，ni为小班i中的样地数，xij为小班i中样地j的观测值，为全研究区样地观测值的均值，为小班i中样地观测值的均值。使用CHM 数据计算时，ni为小班i中的栅格个数，xij为小班i中栅格j的冠层平均高度值，为全研究区冠层平均高的均值，为小班i中冠层平均高的均值。可解释性方差越接近1，则小班内部一致性越高，不同小班间的差异性越大。

3 结果与分析

3.1 定性分析

经过多次试验，对1 m 空间分辨率CHM 数据进行分割时，最优尺度参数为100，形状参数为0.1，紧致度参数为0.5，得到过分割结果。接着对平滑后的5 m 空间分辨率CHM 数据进行分割，最优尺度参数为37，形状参数为0.1，紧致度参数为0.5，得到欠分割结果。

最终区划结果与过分割结果、欠分割结果和人工区划小班的对比如图3 所示。图3-a 中显示过分割结果划分出了一些相对于小班尺度不必要分割开的区域。图3-b 显示欠分割结果没有分割出相对于小班来说不必要的区域，但其边界与小班边缘的符合性较差。对比图3-a 和3-b 可以看出欠分割结果对于小班区域的划分效果比过分割结果好，但过分割结果与小班边缘的符合性更好。优化后的最终区划结果比欠分割结果与林分边缘更相符，比过分割结果更好地区分了不同的林分区域。最终区划结果与人工区划小班的对照如图3-c 所示，结果显示最终区划结果与人工区划小班大体上近似，且最终区划结果与林分边界真实树冠边缘的一致性更高。

图3 基于CHM 的最终区划结果与欠分割结果、过分割结果、人工区划小班对比Fig.3 Under-segmentation results,over-segmentation results,manual sub-compartments and automatic sub-compartments on CHM

3.2 定量分析

3.2.1 形状因子一致性 自动区划结果与参考小班边界一致性的对比如表1 所示。表中各UMA特征从形状和面积等多个角度反映了自动区划结果边界的准确程度，自动区划结果与基于DOM 勾绘的小班、主伐作业小班和人工区划的小班的各特征都较接近，相近程度依次递增。当IoU＞0.7时，自动区划结果与人工区划小班相符性非常好，IoU＞0.5时，自动区划结果与人工区划结果相符性较好。

表1 自动区划结果与参考小班的UMA 特征均值及lOU＞0.7、lOU＞0.5 的自动区划小班占总数的百分比Table 1 The mean values of UMA features and the proportion of automatic sub-compartments with lOU＞0.7、lOU＞0.5compared with reference sub-compartments

由此可知自动区划结果中46%、37%、43%与人工区划小班、主伐作业小班、基于DOM 勾绘的小班相符性非常好，61%、54%、55%与人工区划小班、主伐作业小班、基于DOM 勾绘的小班相符性较好。

交并比较低的小班主要对应以下3 种情况：①人工区划小班边界与地面真实小班边界不符，而自动区划结果边界与真实边界相符。②人工区划小班人为地将同质的林分划分为多个小班，而自动区划时对这部分林分划分的界线与人工区划的不一致。③幼苗幼树小班的树高过低，接近耕地，CHM自动区划方法没有划分开。

3.2.2 林分内部的均质性和外部异质性 自动区划结果的平均胸径可解释性方差为0.97，平均树高可解释性方差为0.98，和人工区划小班相同。可解释性方差精度较高，说明自动区划的小班和人工区划的小班都具有内部一致性高，且不同小班之间区分度好的特点。自动区划结果中小班与样地个数的对应情况与人工区划小班近似，并且二者的可解释性方差一致，说明自动区划方法产生了近似人工区划结果的合理的小班。人工区划小班的平均冠层高度的可解释性方差为83.04%，自动区划小班的平均冠层高度的可解释性方差为84.81%。自动区划小班比人工区划小班的可解释性方差提高了1.77%。说明自动区划小班对不同林分的区分效果比人工区划小班更优。

4 讨论

本研究得到的全场小班区划结果优于人工区划结果，说明CHM 应用于小班区划能够提高边界区划精度，减少主观性。高光谱图像树种分类结果图使小班自动区划结果包含了树种属性。此外，自动区划结果的可解释性方差较高，说明自动区划结果符合小班内部特征基本一致，与相邻小班有明显区别的要求，而定性分析也显示自动区划的小班大部分符合实际情况。本研究方法在约15 000 hm2林场范围的应用时间约为10 h，而人工区划约1 周时间，工作效率显著提高。

在交并比较低的3 种原因中，前2 种是由于人工区划小班结果与地表真实小班的边界不符造成的。这些交并比低的小班，不代表自动区划不合适，而是表现出了自动区划与真实值更加一致，不受主观及其他因素影响的优势。第3 种原因中，在人工区划目视判读时也会有部分幼苗幼树和耕地难以区分开，这些区域需要进行实地调查才能确定。

不同于Dechesne 等[15]的研究需要计算大量的林分特征，本研究的方法可以直接利用高精度的CHM 和树种图进行自动区划来得到和人工小班接近的区划结果。Mustonen 等[13]计算的平均胸径可解释性方差为74%，平均树高可解释性方差为83%。Pukkala 等[16-18]计算的平均树高、平均胸径、断面积可解释性方差为66%～87%，最大树高可解释性方差为84.7%～94.2%。本研究自动区划的小班的平均胸径和平均树高可解释性方差精度相对于这些研究有所提高。本研究的方法在一个完整林场范围内取得了内部一致性精度高且和人工区划小班接近的结果。另外，其他相关区划研究中没有与人工区划小班边界对比，有的林分区划结果虽然绝对精度较高，但图斑较破碎，面积过小，与人工小班差异较大。本研究引入UMA 准则和交并比分析方法，定量计算了使用遥感数据自动区划的小班与人工小班的相符性，并分析了其不同的原因，以及自动区划小班的优势和不足，有利于后续对自动分割小班更好的符合甚至优于人工区划小班的进一步研究。

5 结论

本研究采用机载激光雷达数据提取的CHM 和高光谱影像产生的树种类别图进行了小班自动区划的研究，主要结论如下：（1）本方法充分利用了树高、林分密度、郁闭度等林分结构信息和高精度的树种信息，区划效率远高于人工小班区划；（2）自动区划小班的内部一致性高，外部差异性大，在形状、面积、交并比等多个指标上与人工区划小班相符，自动区划小班的边界较人工区划小班边界更为精准；（3）本方法不需要计算大量特征，适用于林场及更大尺度的小班自动区划。此外，采用机载遥感数据进行自动区划能更好地应用于高精度的森林参数提取等工作，符合小班区划精细化和森林质量精准提升的发展趋势。