CN102419818B - 结合形态学冠层控制和分水岭的LiDAR数据单木提取方法 - Google Patents

Abstract

本发明的目的在于针对CHM、从LiDAR数据中提取单木时，减少对CHM质量的依赖，提高单木提取的精度。具体通过如下方法实现：基于形态学冠层控制，在CHM中确定外部标记；对所述CHM进行全局中值滤波处理，结合所述外部标记，使用局部极值的方式获得中值滤波处理后的CHM的内部标记；结合所述内部标记，对CHM进行两次分水岭操作；利用所述外部标记，分别对分水岭操作结果进行优化；根据单木分割结果，提取单木参数。本发明提供的方法，分割的结果只存在可忽略的过分割现象，使得单木分割的精度得到了提高，且对CHM质量的依赖程度得到降低，具有很强的实用价值。

Description

结合形态学冠层控制和分水岭的LiDAR数据单木提取方法

技术领域

本发明涉及激光雷达数据处理技术，尤其涉及一种结合形态学冠层控制和分水岭的LiDAR数据单木提取方法。 

背景技术

激光雷达(Light Detecting and Ranging，LiDAR)技术是近20年来遥感领域最具革命性的成就之一。根据其特点，LiDAR技术在林业上有着极为广泛的应用，尤其在构建冠层高度模型(Canopy Height Model，CHM)等方面有着突出的优势。 

目前从LiDAR数据中提取单木的算法有很多，其中一些算法是针对CHM进行的。CHM是林业研究中的一个重要模型，它表示了地面以上的冠层高度，反映了植被的垂直高度变化以及植被上表面的水平分布情况。针对CHM、从LiDAR数据中提取单木是指，根据LiDAR数据建立CHM，从CHM中提取单木的方法。目前，从CHM中提取单木的方法主要是区域生长法和标记分水岭法。 

但是，现有的从CHM中提取单木的方法对CHM质量要求较高，且无法保证单木提取的精度。 

发明内容

本发明的目的是提供一种从结合形态学冠层控制和分水岭的LiDAR数据单木提取方法，从而解决在提取单木时对CHM质量的依赖程度过高，以及提取单木精度无法保证的问题。 

本发明的目的通过以下技术方案来实现。 

一种结合形态学冠层控制和分水岭的LiDAR数据单木提取方法，包括： 

步骤1、基于LiDAR数据建立CHM； 

步骤2、基于形态学冠层控制，在所述CHM中确定外部标记； 

步骤3、对所述CHM进行全局中值滤波处理； 

步骤4、结合所述外部标记，使用局部极值的方式获得全局中值滤波处理后的CHM的内部标记； 

步骤5、结合所述内部标记，对步骤1中所述的CHM进行第一次分水岭操作； 

步骤6、利用所述外部标记，剔除第一次分水岭操作结果中的非冠层水槽，得到优化的第一次分水岭操作结果； 

步骤7、查找所述优化的第一次分水岭操作结果中没有内部标记的水槽，对所述没有内部标记的水槽进行重建，得到重建后的CHM； 

步骤8、结合所述内部标记，对重建后的CHM进行第二次分水岭操作； 

步骤9、利用所述外部标记，剔除第二次分水岭操作结果中的非冠层水槽，得到优化的第二次分水岭操作结果； 

步骤10、使用预定的结构元素，对第二次分水岭操作结果进行开操作，得到单木分割结果； 

步骤11、根据单木分割结果，提取单木参数。 

本发明的关键在于始终将冠层控制的概念与单木分割过程结合，并且通过重建，进行了两次的分水岭操作，从而大大的减少了过分割现象。分水岭分割之后，该发明还使用基于自适应结构元素形态学操作对结果进行优化，使得分割的单木冠层更贴近于真实的冠层形状。因此，本发明提供的方法，分割的结果只存在可忽略的过分割现象，使得单木分割的精度得到了提高，且降低了对CHM质量的依赖，具有很强的实用价值。 

附图说明

图1本发明实施例提供的一种方法流程图； 

图2为本发明实施例提供的冠幅测量示意图。 

具体实施方式

本发明提供一种结合形态学冠层控制和分水岭的LiDAR数据单木提取方法，其实现方式如图1所示，包括如下操作： 

步骤1、基于LiDAR数据建立CHM； 

步骤2、基于形态学冠层控制，在步骤1中建立的CHM中确定外部标记； 

所述的外部标记是一个二值图像，用来表示树冠的分布情况。 

所述的外部标记并不参与后续分水岭运算，以最大程度的保留分割的紧凑性，避免为控制过分割而带来的欠分割现象。

外部标记的作用是针对分水岭运算的结果进行冠层控制，以剔除分水岭运算结果中的非冠层水槽。 

可以使用基于形态学的冠层控制进行外部标记的确定。作为举例而非限定，步骤2的具体实现方式可以是： 

步骤21、对上述的CHM进行形态学开操作，开操作的数学表示如下： 

CHM_Open＝CHMοS           公式(1) 

其中，CHM_Open为经过开操作后的CHM，CHM为步骤1中建立的CHM，ο是形态学开操作算子，S为近圆形结构元素。 

步骤22、将开操作的结果与设定的冠层阈值进行比较，根据比较结果获得外部标记的初始二值图像； 

这里需要设定一个冠层阈值threshold_Open，开操作的结果CHM_Open中大于等于threshold_Open的像素被认定是有冠层区域(也就是树冠部分)，而小于该threshold_Open的像素被认定为非冠层区域，即为地表。则步骤22的数学表示如下： 

$\mathrm{outMark}=\left\{\begin{array}{c}1,{\mathrm{CHM}}_{\mathrm{Open}}\&GreaterEqual;{\mathrm{thershold}}_{\mathrm{Open}}\\ 0,{\mathrm{CHM}}_{\mathrm{Open}}<{\mathrm{threshold}}_{\mathrm{Open}}\end{array}\right.$ 公式(2) 

其中，outMark是外部标记的初始二值图像，outMark为1表示冠层区域，outMark为0表示非冠层区域。 

步骤23、通过去除所述外部标记的初始二值图像中小于阈值的独立冠层区域和非冠层区域，得到外部标记二值图像，该外部标记二值图像即外部标记的表现形式。具体可以通过填补法和去除小像素块法实现。 

现有技术中，存在着通过形态学梯度来去除树冠的细节。但CHM上存在较为严重的冠层数据缺陷和丢失，这是本领域技术人员公知的。基于形态学梯度去除树冠细节只能单纯的反应冠层的边界，但无法实现对冠层数据的修复和还原。而本发明提供的方法，基于形态学的冠层控制来确定外部标记(外部标记反映了冠层数据)，然后通过填补法和去除小像素块法对外部标记进行优化，能够保留更多的冠层细节，最大程度的修复和还原冠层的真实形状。 

步骤3、对步骤1中建立的CHM进行全局中值滤波处理，全局中值滤波处理的数学表达如下： 

CHM_Median=Median(CHM)              公式(3) 

其中，CHM_Median是全局中值滤波后的CHM，Median()是中值滤波函数。 

步骤4、结合上述的外部标记，使用局部极值的方式获得全局中值滤波处理后的CHM的内部标记； 

步骤4的数学表示如下： 

$\mathrm{inMark}=\left\{\begin{array}{c}1,\mathrm{CH}{M}_{\mathrm{Median}}(x,y)\&GreaterEqual;\mathrm{Max}(\mathrm{Sub}\_\mathrm{Window}(x,y)),\mathrm{outMark}(x,y)=1\\ 0,{\mathrm{CHM}}_{\mathrm{Median}}(x,y)<\mathrm{Max}(\mathrm{Sub}\_\mathrm{Window}(x,y))\end{array}\right.$ 公式(4) 

其中，inMark为内部标记的二值图像、是内部标记的表现形式，Max()为取最大值函数，Sub_Window()为取(x，y)像素的7*7邻域、不包括中心像素。 

内部标记表示可能的单木位置，相当于表示可能的树顶标记。应当指出的是，单木位置存在于树冠部分，因此，在计算单木位置时，关心的是外部标记取值为1的部分。 

现有技术中，通常假设树顶位于冠层的中间位置。但实际上，树顶往往并不处于冠层的中间位置。发明人在实现本发明的过程中，发现CHM中的单木树顶一定处于一些局部最大值的位置上，但也有可能会出现一个树冠有多个局部最大值的情况。为了确定树顶标记，本发明提供的方法是，首先对CHM进行中值滤波处理，来平滑掉一些同一树冠上的多个树顶的情况；然后结合外部标记和局部极值，查找可能的树顶标记。 

步骤5、结合上述的内部标记，对步骤1中建立的CHM进行第一次分水岭操作； 

在进行第一次分水岭操作之前，需要对步骤1中建立的CHM进行平滑操作，这样可以最大程度的降低分水岭的过分割现象。具体可以使用形态学的闭操作进行平滑处理，其数学表示如下： 

CHM_Close＝CHM●S                    (5) 

其中，CHM_Close是经过形态学闭操作后的CHM，●是形态学闭算子，S为近圆形结构元素。 

步骤5的具体实现方式可以是： 

步骤51、对平滑处理后的CHM进行值反转操作； 

如果基于形态学的闭操作进行平滑处理，则步骤51的具体实现方式是，对CHM_Close进行值反转操作，得到CHM_{Close_Reverse}，其数学表达如下： 

CHM_{Close_Reverse}＝Max(CHM_Close)-CHM_Close+a            公式(6) 

其中，a为一个常数。 

步骤52、在值反转操作后的CHM中加入内部标记，进行第一次分水岭操作。 

步骤6、利用所述的外部标记，剔除第一次分水岭操作结果中的非冠层水槽，得到优化的第一次分水岭操作结果； 

第一次分水岭操作及优化的数学表示如下： 

CHM_Watershed＝Mask(WaterShed(CHM_{Close_Reverse}，inMark)，outMark)    公式(7) 

其中，WaterShed()为分水岭分割函数；Mask()表示使用外部标记的二值图像进行掩膜处理，通过掩膜处理剔除非冠层水槽；CHM_Watershed为优化后的分水岭操作的结果，其表现为许多的水槽。所谓的剔除非冠层水槽是指：将外部标记作为掩膜，如果一个水槽在外部标记的范围之外，则该水槽为非冠层水槽，删除该水槽；如果一个水槽中有部分在外部标记的范围之外，则删除这部分，保留该水槽在外部标记的范围之内的部分。 

步骤7、对优化后的第一次分水岭操作结果进行形态学重建； 

优化后的第一次分水岭操作结果仍然存在过分割现象，因此需要对没有内部标记的水槽进行重建。重建的目的是在第二次分割时过分割水槽将不再被分为一个单独的水槽。 

步骤7的具体实现方式是：查找优化的第一次分水岭操作结果中没有内部标记的水槽，对没有内部标记的水槽进行重建，得到重建后的CHM。以上述数学表达为例，步骤7的具体实现方式可以是：为对CHM_Watershed中的每个水槽进行查找内部标记操作；对于没有内部标记的水槽，则需要将CHM_{Close_Reverse}中该水槽所在像素内的像素值进行重建。重建有多种方式，优选的是使用CHM_{Close_Reverse}中对应水槽的分水岭所在像素的平均值重建水槽像素，其数学表示如下： 

CHM_{Close_Reverse}(Index_{Basins_Without_inMark})＝ 

                                                                    公式(8) 

Mean(CHM_{Close_Reverse}(Boundary_Index_{Basins_Without_inMark})) 

其中，Index_{Basins_Without_inMark}为值没有内部标记的水槽像素索引，Mean()为取均值函数，Boundary_Index_{Basins_Without_inMark}为没有内部标记的水槽边界像素索引。 

步骤8、结合所述内部标记，对重建后的CHM进行第二次分水岭操作； 

步骤9、利用所述外部标记，剔除第二次分水岭操作结果中的非冠层水槽，得到优化的第二次分水岭操作结果； 

优化的第二次分水岭操作结果的分割效果较优化的第一次分水岭操作结果的分割效果有极大的改善，主要体现在绝大部分水槽将会包含内部标记。 

第二次分水岭操作及优化的数学表示如下： 

CHM_Watershed2＝Mask(WaterShed(CHM_{Close_Reverse}，inMark)，outMark)    公式(9) 

其中，公式(9)中的CHM_{Close_Reverse}是经过重建后的CHM_{Close_Reverse}。 

步骤10、使用预定的结构元素，对第二次分水岭操作结果进行开操作，得到单木分割结果； 

经过步骤9得到的优化的第二次分水岭操作结果中存在很多极度不规则的形状，因此需要对水槽进行优化。因为树冠从垂直角度观察是近似圆形，所以可以使用近似圆形的结构元素进行开操作对每个水槽进行优化。结构元素的大小根据水槽的像素数量进行自适应变化。得到的结果将是较接近真实树冠形状的水槽。 

由于绝大部分的树林区域均有严重的重叠和遮挡，造成CHM上冠层形状各异。如果单纯使用步骤9得到的优化的第二次分水岭操作结果作为单木的冠层，结果不够准确。本发明通过步骤10，对每个水槽进行优化，最大程度的还原树冠的真实形状。 

步骤11、根据单木分割结果，提取单木参数。 

通过步骤10得到的单木分割结果，每一个水槽代表的是一棵独立单木的冠层，因此，单木的冠幅可以直接通过水槽区域的计算获得，如图2和公式10所示；结合点云及CHM，可以获得单木的地理位置和树高。 

$\left\{\begin{array}{c}{\mathrm{Radius}}_{\mathrm{average}}=\frac{d_1+d_2+d_3+\&CenterDot;\&CenterDot;\&CenterDot;+d_{n-2}+d_{n-1}+d_n}{n}\\ {\mathrm{Radius}}_{\min }=\mathrm{MIN}(d_1,d_2,d_3,...,d_{n-2},d_{n-1},d_n)\\ {\mathrm{Radius}}_{\max }=\mathrm{MAX}(d_1,d_2,d_3,...,d_{n-2},d_{n-1},d_n)\end{array}\right.$ 公式(10) 

其中，Radius_average是树冠的平均半径，而Radius_min和Radius_max分别是树冠的最小和最大半径，d_n则是图2中的树顶位置到所有边界像素的距离。 

通过上述操作，获得最终的单木分割结果。作为举例而非限定，可以将单木结果输出为矢量文件(*.shp)以供更方便的使用。 

本发明的关键在于始终将冠层控制的概念与单木分割过程，并且通过重建，进行了两次的分水岭操作，从而大大的减少了过分割现象。分水岭分割之后，该发明还使用基于自适应结构元素形态学操作对结果进行优化，使得分割的单木冠层更贴近于真实的冠层形状。因此，本发明提供的方法，分割的结果只存在可忽略的过分割现象，使得单木分割的精度得到了提高，且降低了对CHM质量的依赖程度，具有很强的实用价值。 

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。 

Claims (1)

1.一种结合形态学冠层控制和分水岭的LiDAR数据单木提取方法，其特征在于，包括：

步骤1、基于LiDAR数据建立CHM；

步骤2、对所述CHM进行形态学开操作；将开操作的结果与设定的冠层阈值进行比较，根据比较结果获得外部标记的初始二值图像；通过去除所述外部标记的初始二值图像中小于阈值的独立冠层区域和非冠层区域，得到外部标记二值图像，外部标记用所述外部标记二值图像表示；

步骤3、对所述CHM进行全局中值滤波处理；

步骤4、结合所述外部标记，使用局部极值的方式获得全局中值滤波处理后的CHM的内部标记；

步骤5、对步骤1中所述的CHM进行平滑处理；对平滑处理后的CHM进行值反转操作；在值反转操作后的CHM中加入内部标记，进行第一次分水岭操作；

步骤6、利用所述外部标记，剔除第一次分水岭操作结果中的非冠层水槽，得到优化的第一次分水岭操作结果；

步骤7、查找所述优化的第一次分水岭操作结果中没有内部标记的水槽，对所述没有内部标记的水槽进行重建，得到重建后的CHM；

步骤8、结合所述内部标记，对重建后的CHM进行第二次分水岭操作；

步骤9、利用所述外部标记，剔除第二次分水岭操作结果中的非冠层水槽，得到优化的第二次分水岭操作结果；

步骤10、使用预定的结构元素，对第二次分水岭操作结果进行开操作，得到单木分割结果；

步骤11、根据单木分割结果，提取单木参数。