CN112200083A - 一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法 - Google Patents

Abstract

本发明提供了一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，涉及遥感数据处理技术领域。本发明读取原始机载多光谱LiDAR多波段独立点云数据，形成原始机载多光谱LiDAR多波段独立点云数据集；然后将原始机载多光谱LiDAR独立点云数据集进行异常数据去除及数据融合，形成具有多波段光谱信息的单一点云数据集；进而对具有多波段光谱信息的单一点云数据集提取反应地物类型差异的多光谱强度特征和高程特征；最后将机载多光谱LiDAR数据的分类特征输入到多元高斯混合模型中实现地物聚类，得到每一个数据点的响应度值，按最大响应度原则确定各数据点所属的类别标签，最终得到点云分割结果。

Description

一种基于多元高斯混合模型的机载多光谱LiDAR数据分割 方法

技术领域

本发明涉及遥感数据处理技术领域，尤其涉及一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法。

背景技术

机载激光雷达(Light Detection And Ranging，LIDAR)点云数据分割是点云数据得以应用的先决条件，而自动化、高精度的点云数据分割又将极大扩展点云数据的应用领域。目前采集的点云数据主要来源于单波段机载LiDAR系统，但由于来自LiDAR的后向散射能量取决于目标材料、目标表面粗糙程度和激光波长，因此单波段LiDAR在区分土地覆盖的能力方面受到限制。现有基于单波段的LiDAR点云数据的分割方法主要可分为：(1)将点云数据转换为多次目标识别和分离，通过设定一系列目标识别规则逐步分离各类地物以完成点云分割。这类方法存在误差传递和累积的问题，并且需要反复尝试或凭经验设置各种参数，过程繁琐，对于复杂地物类型指导性不强。(2)将点云数据内插为强度或高程影像，然后从图像中提取统计特征，然后采用基于像元或者对象的分类方法进行地物分类。这种方法将3D点云数据转换为二维栅格数据造成了信息的损失，并且导致了边界模糊、分割不准确的结果。(3)基于高程纹理的分割方法，首先将LiDAR点云内插为高程影像，然后提取高程纹理特征，通过遥感影像分类方法完成分类。这种方法要求待分割地物的高程信息有明显差异，并且单独依靠高程纹理信息进行分类精度不高，需要强度信息进行辅助分类。(4)融合LiDAR点云数据和多光谱或高光谱遥感影像进行分类，可联合利用点云数据的3D空间信息和遥感影像的多光谱信息完成分类。这种方法虽然取得了较好的分类效果，但将不同数据源的数据统一到同一坐标系仍然存在巨大的困难。机载多光谱LiDAR点云是一种新型数据源，其同时包含多光谱和3D空间信息。因此，已有的单波段LiDAR点云分割算法无法直接使用。

现在关于机载多光谱LiDAR的分割方法刚处于起步阶段，主要是对上述第(1)、(2)种方法的改进与应用，依旧未能解决上述存在的问题。基于模型的聚类算法应用统计学的知识，将数据建模成一个概率生成过程，由于这种方法有严密的推导证明及求解算法(Dempster提出的期望最大(EM)算法)，特别是基于高斯混合模型的算法，使得在各种应用中得到了广泛使用，包括信号处理、语音识别、图像分割、高维数据聚类等领域，但在点云分割中应用较少。为了更好的分割机载多光谱LiDAR数据，本发明用高斯混合模型(GaussianMixture Model，GMM)概率分布刻画点云聚类，提出一种基于多元高斯混合模型的机载多光谱LiDAR分割方法，通过将待分割点云数据特征空间看作若干个高斯概率密度函数的线性组合，进而利用GMM拟合机载多光谱LiDAR数据的统计分布，可以实现比单一高斯分布或聚类中心等更复杂的聚类。

本发明探讨了在各个领域广泛应用的高斯混合模型与机载多光谱LiDAR数据提取的特征相结合的可行性与效率，并且综合利用了多光谱LiDAR数据的光谱信息和三维空间信息，有助于其在城市土地覆盖分类的应用。

发明内容

为解决上述技术问题，本发明提出一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，包括以下步骤：

步骤1：读取原始机载多光谱LiDAR数据的各个波段的独立点云数据集，得到原始机载多光谱LiDAR独立点云数据集；

步骤2：对原始机载多光谱LiDAR独立点云数据集进行点云融合，得到具有多波段光谱信息的单一点云数据集；

步骤2.1：将原始机载多光谱LiDAR独立点云数据集中的异常数据进行去除，得到去除异常的多光谱LiDAR独立点云数据集；

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

步骤2.1.2：确定与真实地物对应的最高和最低高程阈值，分别记为T_max和T_min；

步骤2.1.3：将原始机载多光谱LiDAR独立点云数据集中高程低于T_min和高于T_max的激光点判定为高程异常数据，并予以剔除，得到去除高程异常的多光谱LiDAR独立点云数据集；

步骤2.1.4：统计去除高程异常的多光谱LiDAR独立点云数据中各个激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地物对应的强度阈值T_I，将强度大于阈值T_I的点云剔除，得到去除强度异常的多光谱LiDAR独立点云数据集，从而得到去除异常的多光谱LiDAR独立点云数据集；

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

步骤3：对单一点云数据集进行特征提取，构建分割特征向量；

步骤3.1：利用单一点云数据集中各激光点的多波段光谱信息构成光谱特征向量，记为X_B＝[B₁,B₂,...]；其中，B₁，B₂，…分别对应激光点各个波段的激光反射强度值。

步骤3.2：对单一点云数据集进行粗滤波得到地面点集，并由此生成数字高程模型(Digital Elevation Model，DEM)，将单一点云数据集中各激光点的高程减去与激光点平面位置对应的地面高程(由DEM内插得到)得到地物的相对高度，构成各激光点的归一化高程特征向量X_E；

步骤3.2.1：采用不规则三角网的加密滤波算法对单一点云数据集进行滤波得到地面数据集；

步骤3.2.2：通过地面数据集构建规则格网数字高程模型；

步骤3.2.3：将单一点云数据集中每个激光点的高程信息减去其在DEM垂直投影位置处的高程得到地物的高程，构成各激光点的归一化高程特征向量X_E；

步骤4：将不同类型的特征加以组合，得到特征总向量X，并将其作为机载多光谱LiDAR数据的分类特征总向量输入多元高斯混合模型中实现地物聚类分割；

步骤4.1：假设关于分类特征向量X＝[x⁽¹⁾,...,x⁽ⁱ⁾,...,x^(N)]^T的N个样本的集合是独立的，其中，x⁽ⁱ⁾＝[x₁ ⁽ⁱ⁾,...,x_D ⁽ⁱ⁾]为D维特征向量，构建多元高斯混合模型的概率密度函数p(x)，并得到关于特征向量X的对数似然函数Q；

其中，i为数据点索引，i＝1,...,N；x为第i个数据点的特征向量矢量；D为特征向量的维度；X为N×D矩阵；参数K为高斯分布个数，与分割数目相对应；k＝1,...,K为高斯分布的索引；参数μ_k、Σ_k、π_k分别是第k个高斯分布的均值、协方差和权重系数，并且π_k满足

是第k个高斯分布的概率密度函数；对于D维观测变量x，多元高斯分布的形式为

其中，μ为D维均值向量；Σ为协方差矩阵，且是D阶对称正定矩阵；|Σ|为Σ的行列式；

对数似然函数Q表示如下：

步骤4.2：通过EM迭代算法最大化对数似然函数进行聚类，得到各数据点属于各类别的响应度矩阵；

步骤4.2.1：设置分割数目K；将特征向量X作为数据集，对数据集进行聚类分析得到初始分类结果，包括每一类的均值μ和协方差Σ及每个数据点的类标号，通过类标号统计每一类数据点个数，其与总点数的比值作为各类的权重系数π，作为初始的高斯混合模型参数，并计算对数似然函数Q的初始值Q⁽⁰⁾，并设置循环迭代计数器t＝1；

步骤4.2.1.1：在数据集中随机选择一个数据点，将其作为第一个聚类中心c₁；

步骤4.2.1.2：计算所有数据点到当前已有的最近的聚类中心的距离，记为D(x_i)，同时计算每个数据点被选为下一个聚类中心的概率p(x_i)＝D(x_i)/Σ^N _i＝1D(x_i)；选择概率最大的数据点作为新的聚类中心；

步骤4.2.1.3：重复步骤4.2.1.2，直到选出K个聚类中心；

步骤4.2.1.4：计算所有数据点到每个聚类中心的距离，并按照最小距离准则将每个数据点划分到距离其最近的聚类中心；

步骤4.2.1.5：计算每个类中数据点的均值，得到K个新的聚类中心；

步骤4.2.1.6：重复步骤4.2.1.4和4.2.1.5，直到聚类中心不再改变，或者达到预设最大迭代次数；

步骤4.2.2：执行EM算法的E步，依据当前的高斯混合模型参数，计算各高斯分量对数据点的响应度γ(Z_nk)，即第n个数据点属于第k个高斯分布的后验概率；

步骤4.2.3：执行EM算法的M步，使用当前的响应度γ(Z_nk)更新模型参数；

其中，

步骤4.2.4：计算对数似然函数值Q^(t)，检查对数似然函数的收敛性，即计算当前对数似然函数值与上次迭代的对数似然函数值之差||Q^(t)-Q^(t-1)||小于阈值ε，或者迭代次数大于设定阈值T,若满足收敛的条件，则退出循环，否则令t＝t+1，返回步骤4.2.2，直到收敛；

步骤4.3：对响应度矩阵按最大响应度原则确定各数据点所属的类别标签，得到单一点云数据集的分割结果。

本发明所产生的有益效果在于：

本技术方案提供了一种基于多元高斯混合模型的机载多光谱LiDAR数据三维分割方法，该方法首先对原始机载多光谱LiDAR的多波段独立点云数据集进行去噪和融合，得到具有多波段光谱信息的单一点云数据集；然后，通过分析点云数据的光谱信息和高程信息，提取多光谱特征和相对高程特征；最后，利用多元高斯混合模型，通过EM迭代算法最大化对数似然函数进行聚类，得到各数据点属于各地物类别的后验概率，解决机载多光谱LiDAR的三维分割问题。该方法直接基于激光点进行设计，而没有采用像素或体素等数据结构，保留了LiDAR点云数据的真3D优势，避免了像素等数据结构导致的边界模糊、分割不准确的结果；并且利用GMM拟合机载多光谱LiDAR数据的统计分布，将多光谱信息和高程信息有效的结合起来，可以实现更精准的聚类。同时本发明还具有原理直观、易于实现、分割速度快的特点，将有助于机载多光谱LiDAR数据在土地覆盖分类等方面的应用和发展。

附图说明

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

图2为本发明实施例提供的原始机载多光谱LiDAR点云数据；其中，(a)为C1波段强度示意图，(b)为C2波段强度示意图，(c)为C3波段强度示意图；

图3为本发明实施例提供的多元高斯混合模型聚类算法的流程图；

图4为本发明实施例提供的的分割结果；其中，(a)为K＝4的分割结果示意图，(b)为K＝5的分割结果示意图，(c)为K＝6的分割结果示意图。

具体实施方式

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

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

本发明提出了一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，包括以下步骤：

步骤1：读取原始机载多光谱LiDAR数据的各个波段的独立点云数据集，得到原始机载多光谱LiDAR独立点云数据集；

本实施方式中，采用加拿大某公司的Titan机载多光谱LiDAR系统采集数据中的裁剪区域作为实验区域，以检验提出的方法的有效性和可行性。该系统配备了三个不同波段的独立的主动成像通道，波长分别为1550(C1)、1064(C2)和532nm(C3)。其中，每个波段包括采集激光点的三维空间信息及特定波长下地物的激光反射强度信息。实验区域是典型的居民区，包含植被、建筑物、道路和各种附属设施。各波段点云数据的平均点密度约为3.6点/m²。

本实施方式中，如图2所示，记原始机载多光谱LiDAR多波段独立点云数据P_C1,P_C2,P_C3分别为：

其中，C1、C2、C3分别为三个不同波段的标号；i、j、h分别为各个波段激光点的索引，n₁、n₂、n₃分别为各波段激光点的个数，I代表激光点的反射强度值，(x,y,z)对应激光点的三维坐标值。

步骤2：对原始机载多光谱LiDAR独立点云数据集进行点云融合，得到具有多波段光谱信息的单一点云数据集；

步骤2.1：将原始机载多光谱LiDAR独立点云数据集中的异常数据进行去除，得到去除异常的多光谱LiDAR独立点云数据集；

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

步骤2.1.2：确定与真实地物对应的最高和最低高程阈值，分别记为T_max和T_min；

步骤2.1.3：将原始机载多光谱LiDAR独立点云数据集中高程低于T_min和高于T_max的激光点判定为高程异常数据，并予以剔除，得到去除高程异常的多光谱LiDAR独立点云数据集；

步骤2.1.4：统计去除高程异常的多光谱LiDAR独立点云数据中各个激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地物对应的强度阈值T_I，将强度大于阈值T_I的点云剔除，得到去除强度异常的多光谱LiDAR独立点云数据集，从而得到去除异常的多光谱LiDAR独立点云数据集；

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

其中，i’,j’,h’是去除异常数据集中各个波段激光点的索引，n′₁、n′₂、n′₃是去除异常数据集中各波段激光点的个数；

本实施方式中，最高高程阈值T_max、最低高程阈值T_min和强度阈值I根据原始机载多光谱LiDAR点云数据确定。

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

本实施方式中，由于实验区域属于典型的居民区，选择对植被和裸地反射率较强的C1波段作为基础波段，遍历基础波段中的每个激光点，分别搜索半径为r的球体邻域范围内C2、C3两个波段内的激光点N^C2 _pi’，N^C3 _pi’：

根据点云数据的平均密度，设置r＝1m；这里存在两种情况：(1)如果存在多个最近邻点，首先根据激光点的强度值按升序排序，然后根据点个数的奇偶性按下式分别计算此激光点对应的C2、C3波段的强度值；(2)如果不存在最邻近点，则将该波段的强度值设置为0。

记融合后单一点云数据集为：

其中，k为融合点云数据集中各激光点的索引，N为融合点云数据集中激光点的个数。

步骤3：对单一点云数据集进行特征提取，构建分割特征向量；

步骤3.1：利用单一点云数据集中各激光点的多波段光谱信息构成光谱特征向量，记为X_B＝[B₁,B₂,...]；其中，B₁，B₂，…分别对应激光点各个波段的激光反射强度值。

本实施方式中，为消除量纲的影响，将多光谱特征向量X_B离散化至[0,255]。

步骤3.2：对单一点云数据集进行粗滤波得到地面点集，并由此生成数字高程模型(Digital Elevation Model，DEM)，将单一点云数据集中各激光点的高程减去与激光点平面位置对应的地面高程(由DEM内插得到)得到地物的相对高度，构成各激光点的归一化高程特征向量X_E；

步骤3.2.1：采用不规则三角网的加密滤波算法对单一点云数据集进行滤波得到地面数据集；

本实施方式中，首先选取单一点云数据集中的一些高程最低点为初始点构建初始不规则三角形，若单一点云数据集中的点到最近三角形的距离以及单一点云数据集中的点与最近三角形顶点的连线与该三角形的夹角均小于给定的阈值，则将该单一点云数据集中的点加密到该三角网中，并标记为地面数据，通过迭代加密过程，直到没有新点加入三角网时运算终止。完成从单一点云数据集中分离出地面数据集。

步骤3.2.2：通过地面数据集构建规则格网数字高程模型；

步骤3.2.2：通过地面数据集构建规则格网数字高程模型；通过地面点构建规则格网，即将地面点划分为N_x×N_y网格，并将地面点映射到对应的网格中，其中N_x和N_y分别为X方向和Y方向上的网格个数；

本实施方式中，以2倍平均点间距为格网间距，将地面点对应区域划分为N_x×N_y网格，进而将地面点映射到对应的网格中，并将格网内激光点的高程均值作为格网值。

其中，x_max、x_min和y_max、y_min分别代表滤波得到的地面点集中x和y的最大值和最小值；Δx、Δy为x、y方向上的分辨率，其中，S_xy＝{(x_i,y_i),i＝1,...,N}为单一点云数据集P在XOY平面上的投影所得的二维点集，C(S_xy)是点集S_xy的的凸壳，A(C(S_xy))是凸壳C(S_xy)的面积；

步骤3.2.3：将单一点云数据集中每个激光点的高程信息减去其在DEM垂直投影位置处的高程内插值得到地物的高程，构成各激光点的归一化高程特征向量X_E；

本实施方式中，利用反距离加权法(Inverse Distance Weighting，IDW)得到插值点高程。

本实施方式中，为消除量纲的影响，将高程特征向量X_E离散化至[0,255]。

步骤4：将不同类型的特征加以组合，得到特征总向量X，并将其作为机载多光谱LiDAR数据的分类特征总向量输入多元高斯混合模型中实现地物聚类分割，如图3所示；

步骤4.1：假设关于分类特征向量X＝[x⁽¹⁾,...,x⁽ⁱ⁾,...,x^(N)]^T的N个样本的集合是独立的，其中，x⁽ⁱ⁾＝[x₁ ⁽ⁱ⁾,...,x_D ⁽ⁱ⁾]作为数据点，构建多元高斯混合模型的概率密度函数p(x)，并得到关于特征向量X的对数似然函数Q；

其中，i为数据点索引，i＝1,...,N；N为数据点个数；x为第i个数据点的特征向量矢量；D为特征向量的维度；X为N×D矩阵；参数K为高斯分布个数，与分割数目相对应；k＝1,...,K为高斯分布的索引；参数μ_k、Σ_k、π_k分别是第k个高斯分布的均值、协方差和权重系数，并且π_k满足

是第k个高斯分布的概率密度函数；对于D维观测变量x，多元高斯分布的形式为

其中，μ为D维均值向量；Σ为协方差矩阵，且是D阶对称正定矩阵；|Σ|为Σ的行列式。

对数似然函数Q表示如下：

步骤4.2：通过EM迭代算法最大化对数似然函数进行聚类，得到各数据点属于各类别的响应度矩阵；

步骤4.2.1：设置分割数目K；将特征向量X作为数据集，对数据集进行聚类分析得到初始分类结果，包括每一类的均值μ和协方差Σ及每个数据点的类标号，通过类标号统计每一类数据点个数，其与总点数的比值作为各类的权重系数π，作为初始的高斯混合模型参数，并计算对数似然函数Q的初始值Q⁽⁰⁾，并设置循环迭代计数器t＝1；

步骤4.2.1.1：在数据集中随机选择一个数据点，将其作为第一个聚类中心c₁；

步骤4.2.1.2：计算所有数据点到当前已有的最近的聚类中心的距离，记为D(x_i)，同时计算每个数据点被选为下一个聚类中心的概率p(x_i)＝D(x_i)/Σ^N _i＝1D(x_i)；选择概率最大的数据点作为新的聚类中心；

步骤4.2.1.3：重复步骤4.2.1.2，直到选出K个聚类中心；

步骤4.2.1.4：计算所有数据点到每个聚类中心的距离，并按照最小距离准则将每个数据点划分到距离其最近的聚类中心；

步骤4.2.1.5：计算每个类中数据点的均值，得到K个新的聚类中心；

步骤4.2.1.6：重复步骤4.2.1.4和4.2.1.5，直到聚类中心不再改变，或者达到预设最大迭代次数；

步骤4.2.2：执行EM算法的E步，依据当前的高斯混合模型参数，计算各高斯分量对数据点的响应度γ(Z_nk)，即第n个数据点属于第k个高斯分布的后验概率；

步骤4.2.3：执行EM算法的M步，使用当前的响应度γ(Z_nk)更新模型参数；

其中，

步骤4.2.4：计算对数似然函数值Q^(t)(t为迭代次数)，检查对数似然函数的收敛性，即计算当前对数似然函数值与上次迭代计算的对数似然函数值之差||Q^(t)-Q^(t-1)||小于阈值ε，或者迭代次数大于设定阈值T,若满足收敛的条件，则退出循环，否则令t＝t+1，返回步骤4.2.2，直到收敛。

步骤4.3：对响应度矩阵按最大响应度原则确定各数据点所属的类别标签，得到单一点云数据集的分割结果。

本实施方式中，最大响应度原则：对单一点云数据集中的待分割数据点在响应度矩阵中对应的响应度值进行比较，认为最大响应度值所对应的地物类别即为当前数据点所属地物类别：

Z_n＝arg_n{max{γ(z_nk)}}

其中，Z_n表示第n个数据点所属的地物类别，并用Z＝{Z₁,Z₂,…,Z_N}表示数据点的分割结果。

本实施方式中，为了对点云分割结果进行定量的评估，借助Google高分辨率地图，利用商用软件Terrasoild手动对实验数据进行分割，得到各类地物的标准的参考数据；最后，将实验数据和标准数据进行对比，通过混淆矩阵、用户精度(User’s Accuracy)、生产者精度(Producer’s Accuracy)、总体精度(Overall Accuracy)和Kappa系数等评判指标来评估点云分割精度。

本实施方式中，应用本发明提出的方法对实验数据进行分割，分割结果如图4所示。测试数据包含激光点119596个，其中包含有异常数据。经异常数据剔除后，点云个数减少到119374个。经过本发明方法处理后，通过与实际场景目视解译，当聚类类别数K＝4时，将地物分为建筑物、道路、树木和草地；当K＝5时，将地物类别分为建筑物、道路、高植被、中植被和低植被；当K＝6时，将地物类别分为未分类、建筑物、道路、高植被、中植被和低植被。

表1、2、3为本实施例中，聚类类别数K分别为4,5,6时，应用本发明方法对实验数据进行地物分类，对应的分类结果的精度评定。该表中的数据旨在考查不同混合模型的模型分支数对分类结果的影响，并由此确定最优的聚类类别数。

表1 K＝4点云分割结果精度

表2 K＝5点云分割结果精度

表3 K＝6点云分割结果精度

由表1、2、3可知，K＝4、5、6的Kappa系数分别为0.876、0.835和0.827，总体精度分别为90.88％、87.33％和86.62％。这说明：K＝4对应最大的Kappa系数和总体精度，因此，从Kappa系数和总体精度指标来看，K＝4是最佳聚类类别数。通过对比不同类别地物的用户精度和生产者精度，表明：(1)K值得增大不意味着各个地物类别的用户精度和生产者精度必然提高；(2)由于地物存在“同物异谱”和“同谱异物”现象，随着K值的增加，地物类别不断细分，光谱差异较大的同一地物将被细分，光谱差异接近的类别将被错分在一起；同时，聚类类别数的增加也将使地物类别的不断细分。

通过表1、2、3对比可知：最优的地物分割的总体精度及Kappa系数分别为90.88％和0.876。建筑物、道路、树木和草地的用户精度分别为92.35％、87.72％、85.94％、96.24％，生产者精度分别为93.76％、94.98％、93.13％、85.79％。从而验证了本发明提出的方法的有效性。

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

Claims (4)

1.一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，其特征在于，包括以下步骤：

步骤1：读取原始机载多光谱LiDAR数据的各个波段的独立点云数据集，得到原始机载多光谱LiDAR独立点云数据集；

步骤2：对原始机载多光谱LiDAR独立点云数据集进行点云融合，得到具有多波段光谱信息的单一点云数据集；

步骤2.1：将原始机载多光谱LiDAR独立点云数据集中的异常数据进行去除，得到去除异常的多光谱LiDAR独立点云数据集；

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

步骤3：对单一点云数据集进行特征提取，构建分割特征向量；

步骤3.1：利用单一点云数据集中各激光点的多波段光谱信息构成光谱特征向量，记为X_B＝[B₁,B₂,...]；其中，B₁，B₂，…分别对应激光点各个波段的激光反射强度值；

步骤3.2：对单一点云数据集进行粗滤波得到地面点集，并由此生成数字高程模型(Digital Elevation Model，DEM)，将单一点云数据集中各激光点的高程减去与激光点平面位置对应的地面高程(由DEM内插得到)得到地物的相对高度，构成各激光点的归一化高程特征向量X_E；

步骤4：将不同类型的特征加以组合，得到特征总向量X，并将其作为机载多光谱LiDAR数据的分类特征总向量输入多元高斯混合模型中实现地物聚类分割；

步骤4.1：假设关于分类特征向量X＝[x⁽¹⁾,...,x⁽ⁱ⁾,...,x^(N)]^T的N个样本的集合是独立的，其中，x⁽ⁱ⁾＝[x₁ ⁽ⁱ⁾,...,x_D ⁽ⁱ⁾]为D维特征向量，构建多元高斯混合模型的概率密度函数p(x)，并得到关于特征向量X的对数似然函数Q；

其中，i为数据点索引，i＝1,...,N；x为第i个数据点的特征向量矢量；D为特征向量的维度；X为N×D矩阵；参数K为高斯分布个数，与分割数目相对应；k＝1,...,K为高斯分布的索引；参数μ_k、Σ_k、π_k分别是第k个高斯分布的均值、协方差和权重系数，并且π_k满足

是第k个高斯分布的概率密度函数；对于D维观测变量x，多元高斯分布的形式为

其中，μ为D维均值向量；Σ为协方差矩阵，且是D阶对称正定矩阵；|Σ|为Σ的行列式；

对数似然函数Q表示如下：

步骤4.2：通过EM迭代算法最大化对数似然函数进行聚类，得到各数据点属于各类别的响应度矩阵；

步骤4.3：对响应度矩阵按最大响应度原则确定各数据点所属的类别标签，得到单一点云数据集的分割结果。

2.根据权利要求1所述的一种基于多元高斯混合模型的机载多光谱LiDAR分割方法，其特征在于，所述步骤2.1具体包括如下步骤：

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

步骤2.1.2：确定与真实地物对应的最高和最低高程阈值，分别记为T_max和T_min；

步骤2.1.3：将原始机载多光谱LiDAR独立点云数据集中高程低于T_min和高于T_max的激光点判定为高程异常数据，并予以剔除，得到去除高程异常的多光谱LiDAR独立点云数据集；

步骤2.1.4：统计去除高程异常的多光谱LiDAR独立点云数据中各个激光点的强度值的频次，并以直方图的形式可视化显示统计结果；

步骤2.1.5：确定与真实地物对应的强度阈值T_I，将强度大于阈值T_I的点云剔除，得到去除强度异常的多光谱LiDAR独立点云数据集，从而得到去除异常的多光谱LiDAR独立点云数据集。

3.根据权利要求1所述的一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，其特征在于，所述步骤3.2具体包括如下步骤：

步骤3.2.1：采用不规则三角网的加密滤波算法对单一点云数据集进行滤波得到地面数据集；

步骤3.2.2：通过地面数据集构建规则格网数字高程模型；

步骤3.2.3：将单一点云数据集中每个激光点的高程信息减去其在DEM垂直投影位置处的高程得到地物的高程，构成各激光点的归一化高程特征向量X_E。

4.根据权利要求1所述的一种基于多元高斯混合模型的机载多光谱LiDAR数据分割方法，其特征在于，所述步骤4.2具体包括如下步骤：

步骤4.2.1：设置分割数目K；将特征向量X作为数据集，对数据集进行聚类分析得到初始分类结果，包括每一类的均值μ和协方差Σ及每个数据点的类标号，通过类标号统计每一类数据点个数，其与总点数的比值作为各类的权重系数π，作为初始的高斯混合模型参数，并计算对数似然函数Q的初始值Q⁽⁰⁾，并设置循环迭代计数器t＝1；

步骤4.2.1.1：在数据集中随机选择一个数据点，将其作为第一个聚类中心c₁；

步骤4.2.1.2：计算所有数据点到当前已有的最近的聚类中心的距离，记为D(x_i)，同时计算每个数据点被选为下一个聚类中心的概率p(x_i)＝D(x_i)/Σ^N _i＝1D(x_i)；选择概率最大的数据点作为新的聚类中心；

步骤4.2.1.3：重复步骤4.2.1.2，直到选出K个聚类中心；

步骤4.2.1.4：计算所有数据点到每个聚类中心的距离，并按照最小距离准则将每个数据点划分到距离其最近的聚类中心；

步骤4.2.1.5：计算每个类中数据点的均值，得到K个新的聚类中心；

步骤4.2.1.6：重复步骤4.2.1.4和4.2.1.5，直到聚类中心不再改变，或者达到预设最大迭代次数；

步骤4.2.2：执行EM算法的E步，依据当前的高斯混合模型参数，计算各高斯分量对数据点的响应度γ(Z_nk)，即第n个数据点属于第k个高斯分布的后验概率；

步骤4.2.3：执行EM算法的M步，使用当前的响应度γ(Z_nk)更新模型参数；

其中，

步骤4.2.4：计算对数似然函数值Q^(t)，检查对数似然函数的收敛性，即计算当前对数似然函数值与上次迭代的对数似然函数值之差||Q^(t)-Q^(t-1)||小于阈值ε，或者迭代次数大于设定阈值T,若满足收敛的条件，则退出循环，否则令t＝t+1，返回步骤4.2.2，直到收敛。

Images