基于双尺度体元覆盖密度的TLS点云数据单木识别算法
2018-05-14王长青邢万里邢艳秋李梦颖李贝贝
王长青 邢万里 邢艳秋 李梦颖 李贝贝
摘 要：针对现有的地基激光雷达（Terrestrial Laser Scanning， TLS）点云数据单木识别算法存在抗噪性差的问题，本文提出一种基于双尺度体元覆盖密度的TLS点云数据单木识别算法。首先选择内蒙古根河林场的兴安落叶松天然次生林为研究对象，利用徕卡C10三维激光扫描仪获取单测站点云数据；然后通过计算双尺度体元覆盖密度滤除非树干点；最后通过分析体元水平坐标（x，y）位置处的体元z值序列确定滤波后点云数据中的单木位置。研究结果表明：该算法二次滤波结果的平均噪声比为1.66%；滤波后保留的单木数量是实际单木数量的88.94%；滤波后点云数据的单木识别率98.3%，漏检率1.69%，过检率0.56%，实际点云数据的单木识别精度为87.43%。与已有的单测站点云数据单木识别的研究相比，本文提出的单木识别算法简单、抗噪性强且单木识别精度更高，这对于实现复杂密集林分样地单测站点云数据单木的准确识别具有重要意义。
关键词：地基激光雷达TLS， 双尺度， 体元覆盖密度， 滤波， 单木识别
中图分类号：S771.8 文献标识码：A 文章编号：1006-8023（2018）05-0063-09
Abstract： In order to solve the problem of poor noise immunity in the existing individual tree detection algorithm， a individual tree detection algorithm based on voxel cover density of double scale was proposed in this paper. Firstly， the natural secondary forest of Larix gmelinii in Inner Mongolia Genhe Forest was selected as the research object， and the single-test site cloud data was obtained by using the Leica C10 3D laser scanner. Non-trunk points were filtered out by calculated voxel cover density of double scale， and then the individual tree positions were determined by analyzed the z-value sequence in the voxel horizontal coordinated （x，y）. The experiment results showed that the average noise ratio of the secondary filtering result was 1.66%； the number of trees after filtering accounted for 88.94% of the number of all the trees； the individual tree recognition rate was 98.3% in the filtered point cloud data， the missed rate was 1.69%， the more detected rate was 0.56%， and finally the individual tree identification accuracy was 87.43% in the raw data. Compared with the existing research of individual tree identification in single-scan point cloud data， the individual tree detection algorithm of TLS point cloud data based on voxel cover density of double scale proposed in this paper was strong to noise and the detection accuracy of individual tree was high which is very important to realize the accurate individual tree detection of single-scan point cloud data in complex and dense forest stand.
Keywords： Terrestrial Laser Scanning， double scale； voxel cover density； filter； individual tree detection
0 引言
森林是地球上重要的陆地生态系统，在改善全球气候及生态环境等方面均起着重要的作用[1]，森林结构参数的准确估测对于实现林业资源精细化管理及全球碳儲量定量化估测具有重要的意义[2]。
激光雷达作为一种主动遥感技术，可以通过激光测距获取物体表面的高精度三维点云数据，从而实现森林结构参数估测[3-4]。其中，从数据获取方式而言，机载激光雷达（Airborne Laser Scanning， ALS）是从森林冠层上方扫描获得数据，数据质量受森林冠层遮挡影响较大，因而无法获取林冠表面下树冠及树干的细节信息，因此广泛用于林分尺度森林结构参数估测，如林分平均树高、叶面积指数及生物量[5-6]。而地基激光雷达（Terrestrial Laser Scanning， TLS）是从地面向上扫描获得数据，相对而言数据受林冠的影响相对较小且能获得详细的林冠及树干信息，因而，较适于进行森林单木分割[7-8]，从而获得森林单木的详细结构参数，并为其他遥感数据提供高精度的验证数据[9-10]。
单木识别是TLS进行单木分割及进一步实现单木结构参数估测的前提，目前多种算法已成功用于TLS单木识别研究中。按分割对象的不同，先前所用单木分割算法主要可以分为两类，即：基于二维图像的搜索算法和基于三维点云的搜索算法[11]。其中，基于二维图像搜索的单木识别算法是将高程归一化后的切片点云数据转换成灰度图像，然后通过识别灰度图像中的圆来检测单木位置，所用主要算法有先聚类后圆拟合[7，12-13]以及Hough变换[14-16]，但聚类算法无法识别树干和非树干点，同时圆拟合算法容易将部分非树干点拟合出圆，进而导致单木识别率较低；Hough变换法同样会在非树干点位置识别出多个圆，进而导致单木识别精度降低。因此，基于二维图像搜索单木识别算法相对而言受点云噪声点的影响较大，较难实现单木的精确识别。而基于三维点云数据搜索的单木识别方法同样是通过识别树干点的几何特征进行单木识别，常用的算法为圆柱拟合算法，如：Pfeifer等[17]对每个点及距离最近的k个点拟合出一个切面，将平均拟合精度小于2 cm的点视为树干点，然后对剩余的点进行圆柱拟合来识别单木位置。Liang等[18]利用协方差矩阵的特征值和特征向量定义局部坐标系来识别单木树干点，然后通过对识别的树干点进行圓柱拟合进而实现单木识别，对9个半径10 m圆形样地的单测站点云数据的单木识别率为73%，样地平均林分密度为1 022株/hm2。近几年越来越多的研究采用圆柱拟合的方法进行单木识别[18-20]，与基于二维图像的搜索算法相比，圆柱拟合算法先识别树干点后检测单木位置，具有一定的抗噪性，但在树干点识别过程需要对点云数据进行逐点分析，其中设置阈值的难度较大，导致树干点的识别精度较低，而错误保留的非树干点在圆柱拟合过程会降低单木识别精度，因此，改善单木识别算法的抗噪性对于提高单木识别的精度和效率均具有重要意义。
本研究以兴安落叶松天然次生林为研究对象，在单测站扫描点云数据的基础上提出了一种基于双尺度体元覆盖密度的TLS点云数据单木识别算法，以期提高单木识别算法的抗噪性，进而提高单木识别精度，最终实现为其它大尺度遥感数据提供更加准确的地面验证数据。
1 研究区及研究方法
1.1 研究区概况
本研究在内蒙古根河市（120°12'～122°55'E，50°20'～52°30'N）建立研究区，研究区位置如图1所示，该地区位于大兴安岭北段西坡，呼伦贝尔市北部，是中国纬度最高的区域之一，年平均气温-5.3℃，极端低温-58℃，年封冻期210 d以上，海拔高度多在700～1 300 m，平均海拔1 000 m，气候属寒温带湿润型森林气候，并具有大陆季风性气候的某些特征，森林覆盖率91.7%，主要树种有兴安落叶松（Larix gmelinii）、樟子松（Pinus sylvestris）、阔叶树白桦（Betula platyphylla）和山杨（Populus davidiana）等，其中兴安落叶松占林区面积的86.1%，树种总量占大兴安岭所有树种总量的72%。
1.2 数据获取
本研究布设的方形样地和测站位置如图2所示，大样地边长45 m，小样地边长15 m，采用徕卡C10脉冲式三维激光扫描仪进行样地扫描，仪器参数见表1，数据采集时间为2013年8月10号。本文以单测站点云数据为基础数据，选择第1～5、9站距原点水平距离小于15 m的单测站扫描数据进行处理。本研究使用多种仪器手动测量了实验样地的坐标、平均树高、平均胸径等参数，实测数据见表2。
1.3 单木识别算法
本文提出了基于双尺度体元覆盖密度的TLS点云数据单木识别算法。该算法首先对点云数据进行体元化处理，将其分割成两种尺度的体元；然后计算双尺度体元覆盖密度，并对点云数据进行滤波处理，滤除非树干点；最后通过分析体元的z值序列确定单木位置。
1.3.1 点云数据体元化
将所有的点云数据转换为体元，体元边长为l1，记录落入每个体元内的激光点并对其进行体元处理，体元边长为l2，体元坐标的计算方法如公式（1）所示，
式中：（xmin， ymin， zmin）为所有点的x，y，z坐标最小值；l为体元边长；floor为向下取整，此处需尽量避免一个体元内包含两株单木，考虑到单木位置之间的最小间隔距离，所以本文l1设置为0.5 m，l2设置为l1/10=0.05 m。图3（a）为一单木的点云数据被转换为0.5 m的体元，图3（b）则显示了一个0.5 m体元被完整分割成1 000个0.05 m的体元集合。
1.3.2 基于双尺度体元覆盖密度的滤波处理
本文所提基于双尺度体元覆盖密度的单木识别算法在两个体元尺度对点云数据分别进行两次滤波。本文定义体元水平坐标（x，y）位置处沿z轴方向的体元数为该位置的体元覆盖密度。
初次滤波：首先计算每一个0.5 m体元内部0.05 m体元水平位置（x，y）的体元覆盖密度，并将最大体元覆盖密度p1赋值给该0.5 m体元，保留p1≧h1的0.5 m体元。
二次滤波：计算剩余0.5 m体元的水平体元覆盖密度p2，保留p2≧h2位置处的体元。样地调查数据表明，胸径大于5 cm的落叶松树高普遍高于5 m，考虑到单测站扫描数据中由于遮挡造成的树干点云缺失、树干倾斜和初次滤波滤除部分树冠点等因素，本文设置h2=8。针对剩余的0.5 m体元，则统计每个体元水平坐标（x，y）的z值序列，并对二次滤波结果进行以下两个约束：
（1）单木的倾斜生长可能会导致一株单木处识别出2个位置，因此要求z值序列满足zmin≦10和公式（2），
公式（2）中num为计算满足条件的元素数量，如果某个坐标（x，y）位置的z值序列不满足公式（2），则删除该位置沿z轴方向的所有体元。
（2）如果z值序列某序列满足：[a，b]=max（diff（z）），a>8且b<8，则该序列的体元全部被删除，其中diff为计算一次差分序列，max为计算向量的最大值a和最大值对应位置b，即如果z值序列的一次差分序列最大值大于8，且z值序列中
1.3.3 获取单木位置
针对二次滤波剩余的体元z值序列进行分析，对于距离小于2的坐标，本文通过计算两列体元中点云数据之间的距离来判断相邻两列体元是代表1株树还是2株树，相邻两列体元内点云数据之间的距离计算准则为：①如果两个z值序列没有交叉元素的时候，以z值序列的第3个元素对应体元中的点云中心计算距离；②如果交叉元素个数≧3，且交叉元素最小值小于10，以第3个交叉元素对应体元中的点云中心计算距离；③对于交叉元素个数<3和交叉元素最小值大于10的情况，以第1个交叉元素对应体元中的点云中心计算距离。
本研究通过多次实验发现设置单木间距阈值为d=0.35 m时能够尽量减少漏检和过检单木数量，对于相邻两个坐标对应体元中点云数据间距小于d的情况，保留体元覆盖密度最大或点数最多的位置，利用公式（3）計算剩余的体元中心坐标作为单木的位置，式中l=0.5 m。
1.4 精度检验
本文算法的精度检验方法为：在识别的单木位置处构造圆柱点云数据并与原始数据在Terrasolid软件中共同显示，然后通过进行目视判别的方法识别各测站的单木位置及数量，并与本文算法所识别的单木位置进行比较，进而得到本文算法的单木识别精度。样地中单木的相互遮挡会导致单木点云数据不完整，但并不会造成样地中的某株单木被完全遮挡，因此利用目视判别法能够准确地检验本文算法的单木识别精度。
2 结果与分析
2.1 基于双尺度体元覆盖密度的滤波结果
首先对测站1的点云数据进行滤波处理，初次滤波结果如图4所示，其中图4（a）为原始点云数据，图4（b）为初次滤波剩余的点云数据，图4（b）中所有单木的树干都被保留了下来，不存在错误滤除的树干点云数据，但是上层树冠位置也有少量的点云数据被错误保留了下来，这些被错误保留的树冠点分布较为分散，很有可能使得较矮单木被错误识别，因此必须对初次滤波结果作进一步处理。
在初次滤波的基础上进行二次滤波，二次滤波后剩余体元内的点云数据如图5所示，图5中除了黑色多边形内部的树冠点云数据被错误保留之外，其余的树冠数据都被滤除，图中所有单木都有部分点云数据被保留，但是个别单木二次滤波保留的树干明显少于一次滤波，例如图5中黑色边框内部的树干点，该单木二次滤波保留的树干约为一次滤波保留树干的一半，这主要是因为树干倾斜生长导致单木树干贯穿多列体元，而二次滤波只能保留大于体元覆盖密度阈值的体元，所以倾斜生长的单木部分树干会被滤除。
剩余5站点云数据的滤波结果如图6所示，图6（a1～e1）为原始点云数据，图6（a2～e2）为初次滤波后的点云数据，其中包含少量的树冠点和地面点，图6（a3～e3）为二次滤波后的点云数据，各测站点云数据二次滤波结果中，只有少数几株单木树干正上方位置存在被错误保留的离散树冠点。
将滤波后点云数据中的离散树冠点视为噪声并进行统计，6站数据二次滤波后的噪声点统计结果见表3，表中第2站点云数据滤波后的噪声比最大，最大值为2.35%，第5站点云数据滤波后的噪声比最小，最小值为0.8%，6站数据滤波结果的平均噪声比为1.66%。
本文通过在Terrasolid软件中进行目视判别和利用测量工具手动量测的方法统计各测站单木数量，表4为6个测站初次滤波和二次滤波保留的单木数量的统计结果。
表4中二次滤波保留的单木总数是初次滤波的88.94%，除了第一个样地外，其它的样地中二次滤波后的单木数量都小于初次滤波结果后的单木数量，原因有3个：①单测站点云数据中部分单木遮挡严重，大部分树干缺失，不满足体元覆盖密度≧8的约束条件，如图7（a）所示，图中左侧单木高5.26 m，树干缺失2.11 m，剩余点云数据体元覆盖密度<8。②虽然单木树高较大且体元覆盖密度≧8，但z值序列并不满足二次滤波中zmin≦10和公式（2）的约束条件，如图7（a）所示，图中右侧单木缺失了绝大部分树干，尤其是下层树干。③部分单木树高5 m左右，一次滤波后的树干高度虽然>3.5 m，但是由于单木倾斜生长，导致其对应体元坐标（x，y）位置处的体元覆盖密度<8，如图7（b）所示，图中单木树干长3.94 m，图7（c）中单木树干长5.23 m，二者在二次滤波过程都没有被保留下来。表4中第4列为前两个原因导致的单木漏检，漏检率为5.02%，第5列为第3个原因导致的单木漏检，漏检率为6.03%。
2.2 滤波后体元数据的单木识别结果
针对二次滤波后的体元数据，按照1.3.3中的判别条件获取单木位置，统计结果见表5，单木识别结果如图8所示，图中三种单木位置用不同的符号表示并分别标注了序号。
表5中，第1列为二次滤波后各测站的单木数量，第2列为算法识别的单木数量，6站二次滤波数据中共有单木177株，识别单木175株，漏检3株，过检1株，准确识别单木174株，单木准确识别率98.3%，漏检率0.56%，过检率1.69%。第1、2、3、9站点云数据单木准确识别率为100%，第4站漏检2株单木，第5站漏检1株，过检1株。
漏检单木如图9（a）～（c）所示，本文用圆表示识别的单木。造成单木漏检的原因有两个：①未识别的两株单木与相邻的单木为同株双生，二者之间的距离小于阈值0.35 m，如图9（a） （b）所示，图中相邻的两株单木根部相连，这两株漏检单木的位置如图8（d）所示，图中序号1， 2的漏检单木位于序号10， 22的准确识别单木的左侧相邻位置；②两株单木位于同一列体元数据中，如图9（c）所示，图中两株单木位于一列体元中，因此只能识别处一个单木位置，二者的位置如图8（e）所示，图中序号为1的漏检单木和序号为26的准确识别单木坐标重合。
过检单木如图9（d）所示，图中共2株单木但是识别出了3个位置，这3个位置对应图8（e）中序号8， 9的准确识别单木和序号为1的漏检单木，三者坐标沿y轴方向连续分布，错误识别的是中间的单木位置，像这种两株倾斜生长的单木构成的特殊情况，本文提出的单木识别算法沒能准确区分。
3 结论与讨论
本文所提出的基于双尺度体元覆盖密度的TLS点云数据单木识别算法，以兴安落叶松天然次生林为研究对象，以林地单测站点云数据作为基础数据进行实验。现有结论如下：
（1）本文算法具有良好的抗噪性，两次滤波过程可以有效的滤除非树干点，滤波结果的平均噪声比为1.66%，方便了后续的单木位置提取，尤其对天然林复杂环境的点云数据处理具有重要意义。
（2）本研究所选实验样地为中等尺度，样地林分平均密度较高，林地单测站点云数据单木识别率为87.43%，单木识别率较高，实现了密集林分样地单测站点云数据单木的精确识别。
已有的研究中，单站扫描模式下单木检测率介于40%～80%，其中Liang等[18]的研究中的单测站点云数据的单木识别较高，本研究中实验数据半径为15 m，林分平均密度为1 614株/hm2，5 m以上单木的识别率为87.43%，无论是林分密度还是样地半径都要大于Liang等的研究，且识别的最小单木胸径更小（5 cm），但是单木识别精度却高于Liang等的研究。综上所述，本文所提出的基于双尺度体元覆盖密度的TLS点云数据单木识别算法相较于目前主流的圆柱拟合算法，提高了密集林分样地单测站点云数据的单木识别精度。
本文所提出的基于双尺度体元覆盖密度的TLS点云数据单木识别算法抗噪性强且单木识别精度高，但是本研究只分析了落叶松针叶林，在接下来的研究中会尝试对阔叶林或针阔混交林进行处理并验证本算法的单木识别精度。
【参 考 文 献】
[1]BUIZER M， HUMPHREYS D， JONG W D. Climate change and deforestation： The evolution of an intersecting policy domain[J]. Environmental Science & Policy， 2014， 35（1）：1-11.
[2]NEWNHAM G J， ARMSTON J D， CALDERS K， et al. Terrestrial laser scanning for plot-scale forest measurement[J]. Current Forestry Reports， 2015， 1（4）： 239-251.
[3]ZOLKOS S G， GOETZ S J， DUBAYAH R. A meta-analysis of terrestrial aboveground biomass estimation using lidar remote sensing[J]. Remote Sensing of Environment， 2013， 128（1）： 289-298.
[4]TANG H， BROLLY M， ZHAO F， et al. Deriving and validating leaf area index （LAI） at multiple spatial scales through lidar remote sensing： a case study in Sierra National Forest， CA[J]. Remote Sensing of Environment， 2014， 143（3）： 131-141.
[5]N?SSET E， GOBAKKEN T. Estimation of above-and below-ground biomass across regions of the boreal forest zone using airborne laser[J]. Remote Sensing of Environment， 2008， 112（6）： 3079-3090.
[6]SOLBERG S， BRUNNER A， HANSSEN K H， et al. Mapping LAI in a Norway spruce forest using airborne laser scanning[J]. Remote Sensing of Environment， 2009， 113（11）： 2317-2327.
[7]TAO S， WU F， GUO Q， et al. Segmenting tree crowns from terrestrial and mobile LiDAR data by exploring ecological theories[J]. Isprs Journal of Photogrammetry & Remote Sensing， 2015， 110：66-76.
[8]YANG B， DAI W， DONG Z， et al. Automatic forest mapping at individual tree levels from terrestrial laser scanning point clouds with a hierarchical minimum cut method[J]. Remote Sensing， 2016， 8（5）： 372.
[9]JUNG S E， KWAK D A， PARK T， et al. Estimating crown variables of individual trees using airborne and terrestrial laser scanners[J]. Remote Sensing， 2011， 3（11）： 2346-2363.
[10]VINCENT G， ANTIN C， DAUZAT J， et al. Mapping plant area index of tropical forest by Lidar： calibrating ALS with TLS[C]. Silvilaser， 2015.
[11]LIANG X， KANKARE V， HYYPP? J， et al. Terrestrial laser scanning in forest inventories[J]. Isprs Journal of Photogrammetry and Remote Sensing， 2016， 115： 63-77.
[12]KIR?LY G， BROLLY G. 2007. Tree height estimation methods for terrestrial laser scanning in a forest reserve[J]. International Archives of Photogrammetry， Remote Sensing and Spatial Information Sciences， 2007， 36（3/W52）： 211-215.
[13]BU G， WANG P. Adaptive circle-ellipse fitting method for estimating tree diameter based on single terrestrial laser scanning[J]. Journal of Applied Remote Sensing， 2016， 10（2）： 026040.
[14]LINDBERG E， HOLMGREN J， OLOFSSON K， et al. Estimation of stem attributes using a combination of terrestrial and airborne laser scanning[J]. European Journal of Forest Research， 2012， 131（6）： 1917-1931.
[15]OLOFSSON K， HOLMGREN J， OLSSON H. Tree stem and height measurements using terrestrial laser scanning and the RANSAC algorithm[J]. Remote Sensing， 2014， 6（5）： 4323-4344.
[16]劉鲁霞，庞勇，李增元. 2016. 基于地基激光雷达的亚热带森林单木胸径与树高提取[J]. 林业科学， 52（2）：26-37
LIU L X， PANG Y， LI Z Y. Individual tree DBH and height estimation using terrestrial laser scanning （TLS） in a subtropical forest[J]. Scientia Silvae Sinicae， 2016， 52（2）：26-37.
[17]PFEIFER N， GORTE B， Winterhalder D. Automatic reconstruction of single trees from terrestrial laser scanner data[C]. Proceedings of 20th ISPRS Congress， Istanbul： ISPRS， 2004： 114-119.
[18]LIANG X， LITKEY P， HYYPPA J， et al. Automatic stem mapping using single-scan terrestrial laser scanning[J]. IEEE Transactions on Geoscience and Remote Sensing， 2012， 50（2）： 661-670.
[19]LIANG X， KANKARE V， YU X， et al. Automated stem curve measurement using terrestrial laser scanning[J]. IEEE Transactions on Geoscience and Remote Sensing， 2014， 52（3）： 1739-1748.
[20]SRINIVASAN S， POPESCU S C， ERIKSSON M， et al. Terrestrial laser scanning as an effective tool to retrieve tree level height， crown width， and stem diameter[J]. Remote Sensing， 2015， 7（2）： 1877-1896.
