Forest Sub-compartment Delineation Based on Airborne LiDAR Canopy Height Model

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

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

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

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

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

Figure 1.  Canopy height model, tree species map and related data for sub-compartments delineation in Mengjiagang forest farm

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

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

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

Figure 2.  The flowchart of automatic delineation of forest stands based on CHM from airborne LiDAR data and tree species

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

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时滤波结果最优。

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

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

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

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

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

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

其次，采用样地平均胸径、平均树高和5 m × 5 m CHM的冠层平均高数据计算可解释性方差精度，评价小班内部的均质性及不同小班之间的异质性^[16-17,20]。即计算小班内部观测值或冠层平均高的方差之和$ {SS}_{\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{n}} $，各样地值与全研究区均值的方差之和$ {SS}_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}} $的可解释性方差$ {R}^{2} $，公式如下：

使用样地数据计算时，k为小班数，$ {n}_{i} $为小班i中的样地数，$ {x}_{ij} $为小班i中样地j的观测值，$ \stackrel{-}{x} $为全研究区样地观测值的均值，$ \stackrel{-}{{x}_{i}} $为小班i中样地观测值的均值。使用CHM数据计算时，$ {n}_{i} $为小班i中的栅格个数，$ {x}_{ij} $为小班i中栅格j的冠层平均高度值，$ \stackrel{-}{x} $为全研究区冠层平均高的均值，$ \stackrel{-}{{x}_{i}} $为小班i中冠层平均高的均值。可解释性方差越接近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所示，结果显示最终区划结果与人工区划小班大体上近似，且最终区划结果与林分边界真实树冠边缘的一致性更高。

Figure 3.  Under-segmentation results, over-segmentation results, manual sub-compartments and automatic sub-compartments on CHM

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

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

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

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

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

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

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

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