CN111008989A - 一种基于多值体元的机载多光谱lidar三维分割方法 - Google Patents

Abstract

本发明提供一种基于多值体元的机载多光谱LIDAR三维分割方法，涉及遥感数据处理技术领域。本发明步骤如下：步骤1：读取原始机载多光谱LIDAR多波段独立点云数据，形成原始机载多光谱LIDAR多波段独立点云数据集；步骤2：将原始机载多光谱LIDAR多波段独立点云数据集规则化为多值3D体元数据集；步骤3：对多值3D体元数据集进行面向多值3D体元数据分割的3D隐马尔科夫随机场模型模糊聚类。本方法解决了机载多光谱LIDAR数据的点云分割问题，且综合利用了机载多光谱LIDAR数据的多光谱和空间邻域信息，有助于基于多值体元理论的机载LIDAR点云数据处理及应用的发展。

Description

一种基于多值体元的机载多光谱LIDAR三维分割方法

技术领域

本发明涉及遥感数据处理技术领域，尤其涉及一种基于多值体元的机载多光谱LIDAR三维分割方法。

背景技术

点云分割是机载LIDAR数据处理的关键步骤，其是后继的特征描述、目标识别、分类及场景识别等应用的前提。然而目前常用的点云分割方法大多都是针对机载单波段LIDAR点云数据、基于数据的3D空间及反射强度特征设计。机载多光谱LIDAR点云是一种新型数据源，其同时包含多光谱和3D空间信息。因此，已有的单波段LIDAR点云分割算法无法直接使用。若将多光谱LIDAR点云内插成高程、归一化高程或多光谱图像，进而借由传统的图像分割算法实现分割则无法真正发挥多光谱LIDAR点云数据的真3D优势。

按照参与分割的基元类型，现有的基于机载多光谱LIDAR的点云分割方法主要可以划分为：基于点和像元两类方法。基于点基元的分割方法以点或点空间邻域内点集为基元提取反映地物类型差异的几何特征参数，在特征空间通过聚类完成分割。该方法为3D分割算法，但其缺点在于：离散的激光点未明晰表达各点间的邻接和拓扑关系，由此导致分割方法设计困难。为了利用蕴含在点云数据中的空间结构及拓扑信息，采取的以点空间邻域内点集为基元的方式则需研究确定空间邻域尺度的最优值。反复尝试或凭经验设置的解决方案可指导性不强；多尺度空间邻域的解决方案则由于邻域尺度的添加导致特征维数急剧增大，计算复杂度高且耗时。基于像元基元的分割方法将多光谱LiDAR点云内插生成高程或多光谱图像，然后从影像中提取具有高区分度的特征参数并采用基于像素的2D图像分割算法完成分割。该方法将3D点云插值成2D图像进而借由2D图像分割算法完成分割，其优点在于方便算法设计，其缺点在于：一方面存在信息损失、会影响分割结果的准确性，另一方面未真正发挥多光谱LIDAR点云数据的真3D优势。

基于隐马尔科夫随机场模型的模糊聚类(Hidden Markov Random Field ModelBased Fuzzy Clustering,HMRFMBFC)算法可以综合利用多光谱特征和空间邻域信息，但其为2D图像分割算法，无法直接适用于机载多光谱LIDAR数据。

发明内容

本发明要解决的技术问题是针对上述现有技术的不足，提供一种基于多值体元的机载多光谱LIDAR三维分割方法，本方法解决了机载多光谱LIDAR数据的点云分割问题，且综合利用了机载多光谱LIDAR数据的多光谱和空间邻域信息，有助于基于多值体元理论的机载LIDAR点云数据处理及应用的发展。

为解决上述技术问题，本发明所采取的技术方案是：

本发明提供一种基于多值体元的机载多光谱LIDAR三维分割方法，包括以下步骤：

步骤1：读取原始机载多光谱LIDAR多波段独立点云数据，形成原始机载多光谱LIDAR多波段独立点云数据集；

步骤2：将原始机载多光谱LIDAR多波段独立点云数据集规则化为多值3D体元数据集；

步骤2.1：从原始机载多光谱LIDAR多波段独立点云数据中剔除异常数据，得到去除异常的多波段独立点云数据集；

步骤2.2：对去除异常的多波段独立点云数据进行融合，得到具有多波段光谱信息的单一点云数据集；

步骤2.3：将单一点云数据集规则化为多值3D体元数据集；

步骤3：对多值3D体元数据集进行面向多值3D体元数据分割的3D隐马尔科夫随机场模型模糊聚类；

步骤3.1：将多值3D体元数据集中各有值体元作为数据点，构建嵌入三维隐马尔科夫随机场模型模糊聚类方法的目标函数Q；

其中，q为聚类数；m为体元数，u_jk表示第j个体元属于第k类的模糊隶属度，满足0≤u_jk≤1约束；λ为算法的模糊程度；π_jk为约束聚类尺度的变量—先验概率；d_jk为第j个体元属于第k类的非相似性测度；

步骤3.2：利用构建的目标函数对多值3D体元数据进行分割，得到各体元值属于各地物类别的模糊隶属度矩阵；

步骤3.3：对模糊隶属度矩阵按最大隶属度原则进行反模糊化，得到多值3D体元数据的模糊分割结果。

所述步骤2.1具体包括如下步骤：

步骤2.1.1：统计原始机载多光谱LIDAR点云数据中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.2：设置与真实地形对应的最高高程阈值T_h和最低高程阈值T_l；

步骤2.1.3：针对原始机载多光谱LIDAR点云数据中各个激光点，若其高程值高于最高高程阈值T_h或低于最低高程阈值T_l，则该激光点为异常数据，进行剔除，否则保留该激光点，最终获得去除异常的多波段独立点云数据集。

所述步骤2.3具体包括如下步骤：

步骤2.3.1：用单一点云数据集的轴向平行包围盒表示其所覆盖的3D空间范围；

步骤2.3.2：根据单一点云数据集中激光点的平均点间距确定体元在x、y、z方向上的分辨率{△x,△y,△z}，即体元大小；

步骤2.3.3：依据体元分辨率{△x,△y,△z}对轴向平行包围盒进行划分，得到3D体元阵列；

步骤2.3.4：将单一点云数据集中各个激光点映射到3D体元格网，进而根据3D体元格网中包含的激光点的多光谱特征为各体元赋多值，得到多值3D体元数据集；

将含有激光点的体元赋值为激光点的各波段的反射强度值，将不含有激光点的体元赋值为0，所得体元值记做

其中

及

分别表示第j个体元在C₁、C₂、C₃波段的强度值；进一步的，若某一体元中包含多个激光点，则赋值激光点在各波段的反射强度均值；若体元值的各波段反射强度值的强度等级非256，则将赋值后的体元强度值进一步将其离散化到{0,…,255}，然后再进行各体元赋值。

所述步骤3.2具体包括如下步骤：

步骤3.2.1：设置聚类数、模糊程度因子；对多值3D体元数据中有值体元的体元值进行聚类分析得到初始分类结果，将分类结果记入标号场L；利用L求取模糊隶属度矩阵初始值；随机初始化目标函数值Q⁰；

步骤3.2.1.1：将多值3D体元数据集作为样本集，在样本集所在的光谱特征空间随机选取一个体元值作为初始聚类中心；

步骤3.2.1.2：计算各体元值与已有聚类中心的距离，用D(j)表示第j个体元到聚类中心的距离；计算每个体元值被选为下一个聚类中心的概率

选择概率最大的体元值作为新的聚类中心；

步骤3.2.1.3：重复步骤3.2.1.2，直到选出q个聚类中心；

步骤3.2.1.4：按照最小距离准则将各体元值纳入到距离其最近的聚类中心所在的簇，然后，将每个簇的均值作为新的聚类中心点，继续计算每个有值的体元值到聚类中心的距离，重新分类，直到聚类中心不再变化为止；将各体元的分类结果记做标号场L，L＝{l₁,…,l_j,…,l_m}，l_j∈{1,…,k,…,q}是第j个体元的分类标签集；

步骤3.2.2：利用Gibbs分布的隐马尔科夫随机场分布理论，通过定义标号场L的势能函数确定目标函数Q中的先验概率；

步骤3.2.3：求取目标函数中的各体元与聚类中心的非相似性测度；

步骤3.2.4：根据先验概率及非相似性测度，计算各体元属于各地物类别的模糊隶属度，并按最大隶属度原则更新标号场L，即最大隶属度值所对应的地物类别即为当前体元的所属地物类别；

步骤3.2.5：根据模糊隶属度、非相似性测度、先验概率计算目标函数值，记作Q^t，其中t代表迭代次数；

步骤3.2.6：判断是否满足终止条件，若是，则迭代终止，当前模糊隶属度矩阵即为多值3D体元数据中各体元的最优模糊隶属度矩阵；若否，则将当前模糊隶属度作为模糊隶属度初始值，执行步骤3.2.2；

所述终止条件为当前目标函数值与上次计算的目标函数值之差小于阈值ε，或者迭代次数大于设定阈值T；

所述步骤3.2.2具体包括如下步骤：

步骤3.2.2.1：选用Potts模型确定势能函数V_c：对于任一有值体元，比较其分类标签和其空间邻域内体元的分类标签，若相同，则势能为β，否则，势能为0；

步骤3.2.2.2：根据Gibbs分布的势能函数确定目标函数中的先验概率。

先验概率π_jk表示为：

其中，l_j是第j个体元的分类标签，l_jˊ表示其邻域体元的分类标签，N_j为第j个体元的空间邻域体元集合。

所述步骤3.2.3具体包含如下步骤：

步骤3.2.3.1：统计多值3D体元数据集中各聚类的光谱值均值、协方差矩阵；

步骤3.2.3.2：根据所求的光谱值均值、协方差矩阵和先验概率计算各体元与聚类中心的非相似性测度；

第k个聚类的高斯条件概率密度分布函数为p(v_j|l_j＝k)，则用高斯条件概率密度分布函数值的负自然对数定义非相似性测度d_jk；

d_jk＝-logp(v_j|l_j＝k)

其中，w是多光谱特征空间的维数。

采用上述技术方案所产生的有益效果在于：本发明提供的一种基于多值体元的机载多光谱LIDAR三维分割方法，该方法首先将原始机载多光谱LIDAR多波段独立点云数据集规则化为多值3D体元数据集；然后，利用体元值在光谱特征空间呈现的多维正态分布拟合多值3D体元数据中复杂的聚类，并用该后验概率的负对数值表征体元值到聚类中心的非相似性测度。再利用隐马尔可夫随机场模型定义空间邻域体元对中心体元的标号约束，以此为先验概率通过迭代的方式使分割达到“低势能状态(中心体元和空间邻域体元具有相同标号，则势能低)对应的先验概率大”的稳定状态。该方法首先通过规则化的方式将多光谱和3D空间特征融合入同一数据结构，然后将2D图像分割的基于隐马尔科夫随机场模型的模糊聚类经典算法扩展应用至3D，解决了机载多光谱LIDAR数据的点云分割问题，且综合利用了机载多光谱LIDAR数据的多光谱和空间邻域信息，有助于基于多值体元理论的机载LIDAR点云数据处理及应用的发展。

附图说明

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

图2为本发明实施例提供的原始机载多光谱LIDAR点云数据；其中，a为C₁波段点云数据，b为C₂波段点云数据，c为C₃波段点云数据；

图3为本发明实施例提供的将原始机载多光谱LIDAR点云数据规则化为多值3D体元数据集的流程图；

图4为本发明实施例提供的对多值3D体元数据集进行嵌入3D隐马尔科夫随机模型模糊聚类的流程图；

图5为本发明实施例提供的各邻域尺度示意图；其中，a为6邻域，b为18邻域，c为26邻域，d为56邻域，e为80邻域，f为124邻域；

图6为本发明实施例提供的的分割结果。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明以综合利用多光谱和3D空间信息并真正发挥数据的真3D特性为契机，旨在提出面向机载多光谱LIDAR数据的3D点云分割算法，基于规则化后的多值3D体元数据设计实现，创新性地提出了基于多值体元的面向机载多光谱LIDAR点云数据分割的HMRFMBFC 3D算法，将2D图像分割算法中经典的基于隐马尔科夫随机场模型的模糊聚类(Hidden MarkovRandom Field Model Based Fuzzy Clustering,HMRFMBFC)方法扩展至3D，并基于多光谱LIDAR点云规则化后的多值3D体元数据设计实现；

如图1所示，本实施例的方法如下所述。

本发明提供一种基于多值体元的机载多光谱LIDAR三维分割方法，包括以下步骤：

步骤1：读取原始机载多光谱LIDAR多波段独立点云据，形成原始机载多光谱LIDAR多波段独立点云数据集。其中，每个波段包括几何位置及与该位置对应的单波段激光反射强度信息。

本实施方式中，采用某公司采集的机载多光谱LIDAR数据(如图2所示)，以检验提出的方法的有效性和可行性。数据由Titan机载多光谱LIDAR系统获取，传感器在波长分别为1550、1064和532nm等三个成像通道获取LIDAR数据(分别记做C₁、C₂和C₃波段)，观测角度分别为3.5°、0°、7°。数据中包含被道路围绕的住宅建筑物和带有树木附属的居民区。点云数据平均密度约为3.6点/m²。

本实施方式中，记原始机载多光谱LIDAR多波段独立点云数据

为：

其中，C₁、C₂、C₃对应LIDAR点云数据的三个波段；i₁、i₂、i₃是各波段LIDAR点云数据的索引，n₁、n₂、n₃是各波段LIDAR点云数据的个数，I代表反射强度值，p和(x,y,z)分别对应各波段激光点及其对应的坐标。

步骤2：将原始机载多光谱LIDAR多波段独立点云数据集规则化为多值3D体元数据集，具体流程如图3所示。

步骤2.1：从原始机载多光谱LIDAR多波段独立点云数据中剔除异常数据，得到去除异常的多波段独立点云数据集。

步骤2.1.1：统计原始机载多光谱LIDAR点云数据中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果。

步骤2.1.2：确定与真实地形对应的最高高程阈值T_h和最低高程阈值T_l。

步骤2.1.3：针对原始机载多光谱LIDAR点云数据中各个激光点，若其高程值高于最高高程阈值T_h或低于最低高程阈值T_l，则该激光点为异常数据，进行剔除，否则保留该激光点，最终获得去除异常的多波段独立点云数据集。

本实施方式中，去除异常数据集记

如下所示：

其中，i′₁、i′₂、i′₃是去除异常数据集中各波段点云数据的索引，n′、n′₂、n′₃是去除异常数据集中各波段点云数据的个数；

本实施方式中，最高高程阈值T_h和最低高程阈值T_l为常数，其取值需根据原始机载多光谱LIDAR点云数据的空间分布情况确定。本次实验数据中，最高高程阈值T_h＝118，最低高程阈值T_l＝93。

步骤2.2：对去除异常的多波段独立点云数据进行融合，得到具有多波段光谱信息的单一点云数据集。

本实施方式中，对去除异常的多波段独立点云数据中的各个激光点，在半径为r的球体邻域内搜寻其在其它两个波段的空间邻近点，若存在，则依据反距离加权法插值确定该激光点在其它波段的反射强度值，否则，将该激光点在其它波段的反射强度值设置为0。记融合后单一点云数据集为：

其中，i是融合数据中激光点的索引。

本实施方式中，半径r为常数，其取值需根据原始机载多光谱LIDAR点云数据的密度情况确定。

步骤2.3：将单一点云数据集规则化为多值3D体元数据集。

步骤2.3.1：用单一点云数据集的轴向平行包围盒表示其所覆盖的3D空间范围。

本实施方式中，单一点云数据集P的轴向平行包围盒(Aixe Align Bounding Box，AABB)可表示为：AABB＝{(x,y,z)|x_min≤x≤x_max，y_min≤y≤y_max，z_min≤z≤z_max}。其中，(x_max,y_max,z_max)＝max{(x_i,y_i,z_i),i＝1,2,…,n′₁+n′₂+n₃′}，(x_min,y_min,z_min)＝min{(x_i,y_i,z_i),i＝1,2,…,n′₁+n′₂+n′₃}分别代表数据集P中x、y和z坐标的最大值和最小值。

步骤2.3.2：根据单一点云数据集中激光点的平均点间距确定体元在x、y、z方向上的分辨率(△x,△y,△z)，即体元大小。

本实施方式中，体元在x、y、z方向上的分辨率△x、△y、△z的计算公式如下：

其中，S_xy＝{(x_i,y_i),i＝1,...,n₁+n₂+n₃}为单一点云数据集P在XOY平面上的投影所得二维点集，C(S_xy)是点集S_xy的凸壳，A(C(S_xy))是凸壳C(S_xy)的面积。

步骤2.3.3：依据体元分辨率(△x,△y,△z)对轴向平行包围盒进行划分，得到3D体元格网，每一个3D体元格网单元即为体元。

本实施方式中，基于体元分辨率(△x,△y,△z)即可将AABB划分为3D体元格网，用3D体元阵列表示。基于体元分辨率(△x,△y,△z)对单一点云数据集的轴向平行包围盒划分为3D体元格网，用3D体元阵列表示，设V是3D体元阵列中的体元集合：

V＝{v₁、v₂、…、v_j、…、v_m}

其中，m是体元数，j＝1,…,m；v_j＝{r_j,c_j,l_j}是第j个体元的体元值，r_j、c_j、l_j分别代表第j个体元在3D体元阵列中的行坐标、列坐标和层号；X方向上的体元数R、Y方向上的体元数C及Z方向上的体元数L由下式确定：

其中，x_max、y_max、z_max分别代表单一点云数据集中x、y和z坐标的最大值，x_min、y_min、z_min分别代表单一点云数据集中x、y和z坐标的最小值，

为向上取整操作符，则体元数m＝R*C*L；

步骤2.3.4：将单一点云数据集中各个激光点映射到3D体元格网，进而根据3D体元格网中包含的激光点的多光谱特征为各体元赋多值，得到多值3D体元数据集。

本实施方式中，依下式，将单一点云数据集P中各个激光点映射到3D体元格网，进而根据3D体元格网中包含的激光点的多光谱特征为各体元赋值。其中，将含有激光点的体元赋值为激光点的各波段的反射强度值、将不含有激光点的体元赋值为0。

其中，

为向下取整操作符。若某一体元中包含多个激光点，则赋值激光点在各波段的反射强度均值。进一步，若体元值的各波段反射强度值的等级非256，则将其离散化到{0,…,255}，所得体元值记做

由此，得到多值3D体元数据集，完成对单一点云数据集的规则化。

步骤3：对多值3D体元数据集进行面向多值3D体元数据分割的3D隐马尔科夫随机场模型模糊聚类。具体流程如图4所示。

步骤3.1：将多值3D体元数据集中各有值体元作为数据点，构建嵌入三维隐马尔科夫随机场模型模糊聚类方法(即3D HMRFMBFC)的目标函数Q:

其中，q为聚类数；m为体元数，u_jk表示第j个体元属于第k类的模糊隶属度，满足0≤u_jk≤1约束；λ为算法的模糊程度；π_jk为约束聚类尺度的变量—先验概率；d_jk为第j个体元属于第k类的非相似性测度；

本实施方式中，各体元的分类标签的随机性，利用Gibbs分布的马尔科夫随机场(Markov Random Field，MRF)分布理论、通过定义标号场的势能函数确定目标函数Q中的先验概率π_jk。d_jk为第j个体元属于第k类的条件概率，即观测随机场，代表了多值3D体元数据中某一类别体元值的分布情况，通常，HMRF状态的观测变量的条件概率分布是多维高斯正态分布，可利用该概率的负对数值表征体元值到聚类中心的非相似性测度d_jk。该HMRFMBFC目标函数参数包括模糊程度因子、模糊隶属度矩阵、先验概率、非相似性测度等参数。通过定义标号场的势能函数确定目标函数中的先验概率，以考虑体元的空间邻域关系对聚类结果的影响；通过定义非相似性测度考虑光谱特征对聚类结果的影响，因此，该方法同时考虑了3D空间和光谱信息对聚类结果的影响。

步骤3.2：利用构建的目标函数对多值3D体元数据进行分割，得到各体元值属于各地物类别的模糊隶属度矩阵；

步骤3.2.1：设置聚类数(即机载多光谱LIDAR数据中的地物类别个数)、模糊程度因子(即分割结果的混乱程度)；对多值3D体元数据中有值体元的体元值进行聚类分析得到初始分类结果，将分类结果记入标号场L；利用L求取模糊隶属度矩阵初始值；随机初始化目标函数值Q⁰。

本实施方式中，聚类数q取值为4；模糊程度因子λ取值1.5。

本实施方式中，聚类分析算法采用K-means++，具体步骤如下：

步骤3.2.1.1：将多值3D体元数据中的各有值体元的体元值

(即多光谱特征)作为样本集，在样本集所在的光谱特征空间随机选取一个体元值作为初始聚类中心。

步骤3.2.1.2：首先，计算各体元值与已有聚类中心的距离，用D(j)表示第j个体元到聚类中心的距离，然后，计算每个体元值被选为下一个聚类中心的概率

最后，选择概率最大的体元值作为新的聚类中心。

步骤3.2.1.3：重复步骤b直到选出q个聚类中心。

步骤3.2.1.4：首先，按照最小距离准则将各体元值纳入到距离其最近的聚类中心所在的簇，然后，将每个簇的均值作为新的聚类中心点，继续计算每个有值体元值到聚类中心的距离，重新分类直到聚类中心不再变化为止。

将各体元的分类结果记做标号场L。L＝{l₁,…,l_j,…,l_m}，l_j∈{1,…,k,…,q}是第j个体元的分类标签，若分成4类，则l_j∈{1,2,3,4}。

本实施方式中，隶属度矩阵u_jk的生成方案：第j个体元的所属地物类别赋予q类中的最大隶属度值，且满足

步骤3.2.2：利用Gibbs(吉布斯)分布的隐马尔科夫随机场(Hidden MarkowRandom Field，HMRF)分布理论、通过定义标号场L的势能函数确定目标函数Q中的先验概率；

步骤3.2.2.1：选用Potts模型定义势能函数V_c：对于任一有值体元，比较其分类标签和其空间邻域内体元的分类标签，若相同，则势能为β，否则，势能为0。

在本实施方式中，在标号场中定义体元间空间邻域关系、势能函数公式如下：

其中，l_j是第j个体元的分类标签，l_jˊ表示其邻域体元的分类标签，在6、18、26、56、80及124等不同空间邻域(如图5所示)内定义势能函数，当中心体元与邻域体元的分类标签相同时，达到稳定状态，势能为0，否则，势能为1，β为势能函数中的耦合系数，β由经验获得，β∈[0.5,1]。当空间邻域为6时，式中V_c∈{0,1,…,6}，其它空间领域尺度的V_c类推。

本实施方式中，应用不同的空间邻域尺度会得到不同的势能，并进而影响分割结果的精度。最佳邻域尺度将在实验中确定。

步骤3.2.2.2：根据Gibbs分布的势能函数确定目标函数中的先验概率。

本实施方式中，先验概率π_jk由下式确定。

其中，N_j为第j个体元的空间邻域体元集合。

步骤3.2.3：求取目标函数中的各体元与聚类中心的非相似性测度。

步骤3.2.3.1：统计多值3D体元数据中各聚类的光谱值均值、协方差矩阵。

本实施方式中，设μ_k＝(μ_k1,μ_k2,μ_k3)^T为第k个聚类的均值，μ_k1,μ_k2,μ_k3是第k个聚类均值的3个分量；∑_k＝(∑_k1,∑_k2,∑_k3)^T为第k个聚类的协方差矩阵，∑_k1,∑_k2,∑_k3是第k个聚类协方差矩阵的3个分量。

步骤3.2.3.2：根据所求的光谱值均值、协方差矩阵和先验概率计算各体元与聚类中心的非相似性测度。

本实施方式中，记第k个聚类的高斯条件概率密度分布函数为p(v_j|l_j＝k)，则可以用高斯条件概率密度分布函数值的负自然对数定义非相似性测度d_jk。

其中，w是多光谱特征空间的维数，w＝3。

步骤3.2.4：根据先验概率及非相似性测度，计算各体元属于各地物类别的模糊隶属度，并按最大隶属度原更新标号场L(即认为最大隶属度值所对应的地物类别即为当前体元的所属地物类别)。

本实施方式中，当前模糊隶属度记做

由下式确定。

步骤3.2.5：根据模糊隶属度、非相似性测度、先验概率计算HMRFMBFC目标函数值。

本实施方式中，当前目标函数值记做Q¹，由下式确定。

步骤3.2.6：若当前目标函数值与上次计算的目标函数值之差小于阈值ε，或者迭代次数大于设定阈值T，则迭代终止，当前模糊隶属度矩阵即为多值3D体元数据中各体元的最优模糊隶属度矩阵。否则将当前模糊隶属度作为模糊隶属度初始值，回到步骤3.2.2。

在本实施方式中，如果|Q¹-Q⁰|≥ε，则令Q¹＝Q⁰，重复步骤3.3～3.7；若|Q¹-Q⁰|<ε迭代停止，输出当前模糊隶属度

作为最优模糊隶属度矩阵。阈值ε＝0.0005，最大迭代次数T＝100。

步骤3.3：对上述模糊隶属度矩阵按最大隶属度原则进行反模糊化，，得到多值3D体元数据的模糊分割结果，实现多值3D体元数据的模糊分割。

本实施方式中，最大隶属度原则：对多值3D体元数据中的待分割有值体元在模糊隶属度矩阵中对应的模糊隶属度值进行比较，认为最大隶属度值所对应的地物类别即为当前体元的所属地物类别：

Z_j＝arg_j{max{u_jk}}

其中，Z_j表示第j个体元所属的地物类别，并用Z＝{Z₁,Z₂,…,Z_m}表示多值3D体元数据的模糊分割的结果。

本发明可以在CPU AMD A10 PRO-7800B R7 3.50Ghz、内存4GB、Windows 7旗舰版系统上使用MATLAB R2016b软件编程实现仿真。

本实施方式中，采用商用Terrasoild软件、手工将测试数据准确分类为建筑物、树木、草地、道路等4类作为标准参考数据，以定量评价本发明提出的方法的精度。

本发明提出的HMRFMBFC方法的成果是以体元形式表示的，而参考数据中则是以离散的LIDAR激光点云数据表达的。为和参考数据作对比以评价本发明提出的方法的精度，首先统计本方法所分割出的各类体元内包含的原始机载LIDAR点云数据，然后和参考数据进行比对，进而用正确率(正确分类的激光点数占检测结果中激光点总数的比例)、完整度(正确分类的激光点数占标准数据中激光点总数的比例)、总体精度、混淆矩阵及Kappa系数来定量评价本发明所提出的HMRFMBFC方法的有效性。

表1为本实施例中，邻域尺度分别为6、18、26、56、80和124时，应用本发明方法对测试数据进行地物分类，对应的分类结果的Kappa系数。该表中的数据旨在考查不同领域尺度对分类结果的影响，并由此确定最优空间邻域尺度。

表1不同空间邻域尺度的HMRFMBFC算法分割精度

由表1可知，6、18、26、56、80和124邻域的Kappa系数分别为57.0％、60.9％、61.3％、76.0％、83.0％和73.9％。这说明：(1)80邻域对应最大的Kappa系数，因此，从Kappa系数指标来看，80邻域是最佳空间邻域尺度。(2)空间邻域尺度的增加并不意味着检测精度的必然提高。本发明提出的隐马尔科夫随机场模型的思想是，某个体元的分类情况只与其空间邻域体元的分类情况有关，而与其它体元的分类信息无关。即可用空间邻域体元的分类情况预测中心体元属于各类的先验概率。但是，若空间邻域尺度过小，则邻域系统无法完整地反映地物的空间分布信息，预测的低精度必然影响分割精度，如6、18、26等邻域的HMRFMBFC算法的Kappa系数所示，其值较小。随着空间邻域尺度的增大，当邻域系统正确反映周围地物目标点空间分布时，则可实现中心体元的标号的正确预测，从而达到最高的分割精度。如80邻域的HMRFMBFC算法的Kappa系数所示，达到最大值。空间邻域尺度过大时，邻域系统纳入错误的地物空间分布信息的概率增大，预测出错的概率较大，分割精度同样较低。这可以解释为何124邻域对比80邻域的精度反而有所下降。

本实施方式中，应用本发明提出的方法对测试数据进行分割，分割结果如图6所示。测试数据包含激光点118596个，其中包含有异常数据。经异常数据剔除后，点云个数减少到118386个。上述点云数据被规则化为多值3D体元数据(体元分辨率为0.6m×0.6m×0.6m、3D体元数据尺寸为398×400×42)，其中包含83511个非0值体元。经过HMRFMBFC处理，从上述体元中分离出第Ⅰ、Ⅱ、Ⅲ及Ⅳ类的体元分别为18413、20070、31667及13361个。统计各类体元内包含的激光点数，分别为19210、46826、23931及28419个。

表2、表3为本实施例中，应用本发明提出的方法以参考数据为标准对测试数据的80邻域尺度下的点云分割精度进行的定量评价。

表2点云分割结果的精度

实验数据	总体精度(％)	Kappa系数(％)
			Area 1	87.9	83.0

表3 Area1分割结果混淆矩阵

由表2和表3可知：第Ⅰ、Ⅱ、Ⅲ及Ⅳ类地物分割的平均总体精度及Kappa系数分别为87.9％和83.0％。第Ⅰ、Ⅱ、Ⅲ及Ⅳ类地物的完整率分别为91.8％、88.9％、85.3％、88.3％，正确率分别为81.5％、82.4％、93.2％、91.2％。从而验证了本发明提出的方法的有效性。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (6)

1.一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：包括以下步骤：

步骤1：读取原始机载多光谱LIDAR多波段独立点云数据，形成原始机载多光谱LIDAR多波段独立点云数据集；

步骤2：将原始机载多光谱LIDAR多波段独立点云数据集规则化为多值3D体元数据集；

步骤2.1：从原始机载多光谱LIDAR多波段独立点云数据中剔除异常数据，得到去除异常的多波段独立点云数据集；

步骤2.2：对去除异常的多波段独立点云数据进行融合，得到具有多波段光谱信息的单一点云数据集；

步骤2.3：将单一点云数据集规则化为多值3D体元数据集；

步骤3：对多值3D体元数据集进行面向多值3D体元数据分割的3D隐马尔科夫随机场模型模糊聚类；

步骤3.1：将多值3D体元数据集中各有值体元作为数据点，构建嵌入三维隐马尔科夫随机场模型模糊聚类方法的目标函数Q；

其中，q为聚类数；m为体元数，u_jk表示第j个体元属于第k类的模糊隶属度，满足0≤u_jk≤1约束；λ为算法的模糊程度；π_jk为约束聚类尺度的变量—先验概率；d_jk为第j个体元属于第k类的非相似性测度；

步骤3.2：利用构建的目标函数对多值3D体元数据进行分割，得到各体元值属于各地物类别的模糊隶属度矩阵；

步骤3.3：对模糊隶属度矩阵按最大隶属度原则进行反模糊化，得到多值3D体元数据的模糊分割结果。

2.根据权利要求1所述的一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：所述步骤2.1具体包括如下步骤：

步骤2.1.1：统计原始机载多光谱LIDAR点云数据中各个激光点高程值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.2：设置与真实地形对应的最高高程阈值T_h和最低高程阈值T_l；

步骤2.1.3：针对原始机载多光谱LIDAR点云数据中各个激光点，若其高程值高于最高高程阈值T_h或低于最低高程阈值T_l，则该激光点为异常数据，进行剔除，否则保留该激光点，最终获得去除异常的多波段独立点云数据集。

3.根据权利要求1所述的一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：所述步骤2.3具体包括如下步骤：

步骤2.3.1：用单一点云数据集的轴向平行包围盒表示其所覆盖的3D空间范围；

步骤2.3.2：根据单一点云数据集中激光点的平均点间距确定体元在x、y、z方向上的分辨率{△x,△y,△z}，即体元大小；

步骤2.3.3：依据体元分辨率{△x,△y,△z}对轴向平行包围盒进行划分，得到3D体元阵列；

步骤2.3.4：将单一点云数据集中各个激光点映射到3D体元格网，进而根据3D体元格网中包含的激光点的多光谱特征为各体元赋多值，得到多值3D体元数据集；

将含有激光点的体元赋值为激光点的各波段的反射强度值，将不含有激光点的体元赋值为0，所得体元值记做

其中

及

分别表示第j个体元在C₁、C₂、C₃波段的强度值；进一步的，若某一体元中包含多个激光点，则赋值激光点在各波段的反射强度均值；若体元值的各波段反射强度值的强度等级非256，则将赋值后的体元强度值进一步将其离散化到{0,…,255}，然后再进行各体元赋值。

4.根据权利要求1所述的一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：所述步骤3.2具体包括如下步骤：

步骤3.2.1：设置聚类数、模糊程度因子；对多值3D体元数据中有值体元的体元值进行聚类分析得到初始分类结果，将分类结果记入标号场L；利用L求取模糊隶属度矩阵初始值；随机初始化目标函数值Q⁰；

步骤3.2.1.1：将多值3D体元数据集作为样本集，在样本集所在的光谱特征空间随机选取一个体元值作为初始聚类中心；

步骤3.2.1.2：计算各体元值与已有聚类中心的距离，用D(j)表示第j个体元到聚类中心的距离；计算每个体元值被选为下一个聚类中心的概率

选择概率最大的体元值作为新的聚类中心；

步骤3.2.1.3：重复步骤3.2.1.2，直到选出q个聚类中心；

步骤3.2.1.4：按照最小距离准则将各体元值纳入到距离其最近的聚类中心所在的簇，然后，将每个簇的均值作为新的聚类中心点，继续计算每个有值的体元值到聚类中心的距离，重新分类，直到聚类中心不再变化为止；将各体元的分类结果记做标号场L，L＝{l₁,…,l_j,…,l_m}，l_j∈{1,…,k,…,q}是第j个体元的分类标签集；

步骤3.2.2：利用Gibbs分布的隐马尔科夫随机场分布理论，通过定义标号场L的势能函数确定目标函数Q中的先验概率；

步骤3.2.3：求取目标函数中的各体元与聚类中心的非相似性测度；

步骤3.2.4：根据先验概率及非相似性测度，计算各体元属于各地物类别的模糊隶属度，并按最大隶属度原则更新标号场L，即最大隶属度值所对应的地物类别即为当前体元的所属地物类别；

步骤3.2.5：根据模糊隶属度、非相似性测度、先验概率计算目标函数值，记作Q^t，其中t代表迭代次数；

步骤3.2.6：判断是否满足终止条件，若是，则迭代终止，当前模糊隶属度矩阵即为多值3D体元数据中各体元的最优模糊隶属度矩阵；若否，则将当前模糊隶属度作为模糊隶属度初始值，执行步骤3.2.2；

所述终止条件为当前目标函数值与上次计算的目标函数值之差小于阈值ε，或者迭代次数大于设定阈值T。

5.根据权利要求4所述的一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：所述步骤3.2.2具体包括如下步骤：

步骤3.2.2.1：选用Potts模型确定势能函数V_c：对于任一有值体元，比较其分类标签和其空间邻域内体元的分类标签，若相同，则势能为β，否则，势能为0；

步骤3.2.2.2：根据Gibbs分布的势能函数确定目标函数中的先验概率；

先验概率π_jk表示为：

其中，l_j是第j个体元的分类标签，l_jˊ表示其邻域体元的分类标签，N_j为第j个体元的空间邻域体元集合。

6.根据权利要求4所述的一种基于多值体元的机载多光谱LIDAR三维分割方法，其特征在于：所述步骤3.2.3具体包含如下步骤：

步骤3.2.3.1：统计多值3D体元数据集中各聚类的光谱值均值、协方差矩阵；

步骤3.2.3.2：根据所求的光谱值均值、协方差矩阵和先验概率计算各体元与聚类中心的非相似性测度；

第k个聚类的高斯条件概率密度分布函数为p(v_j|l_j＝k)，则用高斯条件概率密度分布函数值的负自然对数定义非相似性测度d_jk；

d_jk＝-log p(v_j|l_j＝k)

其中，w是多光谱特征空间的维数。

Images