基于多级曲率约束的高斯混合模型点云滤波方法与流程

[0001]
本发明属于激光雷达应用、计算机视觉和地理信息系统领域，主要包括机载lidar点云的地面滤波过程中的点云去噪、曲面格网构建以及根据后验概率确定点的类别等方面。


背景技术：

[0002]
机载lidar系统通过测量激光回波的时间来得到目标至仪器中心的距离。只要给定地理参考系，系统能够快速得到高精度的点云数据。与传统由摄影测量得到点云的方法比较， lidar系统的优势在于较短的处理时间以及更高的测量精度。当前，lidar技术已经被推广到多个领域，如数字地形模型构建、三维建筑模型构建、路面提取以及森林站参数的估计等等。在众多的应用中，将地面测量与非地面测量在原始点云中区分开(滤波)是一项必不可少的工作。同样，精确的数字高程模型(dem)构建也需要预先将非地面点移除，进而由地学插值方法计算得到栅格格网。
[0003]
机载lidar技术革新了获取高精度dem的方式，这主要是因为其获取点云效率高、精度高以及价格低廉。相较于数字摄影测量技术，机载lidar提供了一种对大范围区域的高精度地形测绘方案，许多研究者将注意力集中于如何使用lidar点云数据获取地形信息。为了获取dem，必须要进行的步骤是过滤掉非地面点云，包括来自桥梁、建筑、植被等的测量值。值得注意的是，虽然机载lidar技术能够在较短的时间内获得庞大的数据量，但是其滤波以及手动质量控制却要占到整个后处理时间的60％-80％。因此，如何减少点云数据处理算法的时间复杂度并减少其参数调整时间是当前研究工作的主要问题。
[0004]
在过去的二十年间，有数十种地面滤波算法被提出来。然而，从处理时间以及处理精度来看，滤波依旧是点云数据处理的一大挑战。sithole和vosselman在2004年对8种经典的滤波算法于农村以及城市试验区数据开展了综合性的比较实验工作。他们指出所有的算法在较为平坦的地形上都能够表现出色。然而，当地形变得陡峭并伴随有复杂的地物(如桥梁、低矮灌木等)，滤波器将会失效而表现得不尽人意。总体而言，基于曲面的滤波器在大多数场景的表现大多数较好。这主要是因为该类型滤波器相较于其他类型滤波器利用了更多的点云数据上下文信息。自此，他们选用的位于农村及城市的15块测试样本被用作测试地面滤波算法的标准数据集。
[0005]
机载lidar点云的地面滤波算法主要可以分为四大类：基于坡度的滤波器，基于数学形态学的滤波器，基于曲面插值的滤波器，基于分割策略的滤波器。其中基于数学形态学以及基于曲面插值的滤波器是广泛使用的方法。特别地，商用软件terrascan集成了经典渐进三角加密滤波(ptd)算法，能够以相对高的精度得到地形测量。但是，该算法的一大不足就是需要复杂的参数调整，且因三角面片会在特殊地形上变形严重而无法满足任意地形的地面滤波精度要求。而基于坡度的滤波器多假设地形点测量间的坡度将小于非地形点，通过给定坡度阈值滤除非地形点，由于假设条件大多适用于较为平坦的地形，而大多数的地形情况是复杂多变的，因此限制了该类方法在实际工作中的应用。基于分割策略的滤波器
的基本思想是先聚类后分类，通过同质性原则采用聚类算法如区域增长方法进行测量点的合并，然后给定角度或者距离阈值对各个分割段分类为地面段与非地面段，从而达到滤波目的，但该方法通常需要选择合适的聚类方法保证分割段内部均质性，且伴随着复杂参数设置问题，在滤波工作中较少使用。
[0006]
现今的大多数滤波算法都需要根据具体的试验区域进行复杂的参数调整才能达到要求的精度。然而，这对于许多测绘业界工作者而言这并不是十分友好，因为他们需要对算法机制有着一定的了解且需要大量的实验调整参数。举例而言，当使用基于数学形态学的滤波器时，算法往往假定坡度s为一恒定的常数。但这个假设在实际工作中无法得到满足，因为受到地表内、外力作用下，地形坡度总是不断变化的。使用恒定的s值可能会导致，高度阈值dt计算错误，最终产生大量的遗漏误差或者委托误差，这在后续dem构建中将导致灾难性错误。


技术实现要素：

[0007]
本发明的目的在于针对现有地面滤波算法存在的复杂参数调整、滤波结果需要后期交互验证修正的不足，提出一种新型的基于多级曲率约束的高斯混合模型点云滤波方法。
[0008]
基于多级曲率约束的高斯混合模型点云滤波方法，包括以下步骤：
[0009]
步骤一，输入机载lidar点云，根据大数定律，首先计算每个点云到其k近邻点的平均距离，然后基于大于m倍标准差为粗差的原则，剔除不满足条件的点云，得到待处理点云；
[0010]
步骤二，对待处理点云所在区域建立网格索引，并由局部最小值方法获得的地面种子点，利用薄板样条函数(tps)进行插值，剔除与插值曲面高差大于曲率阈值的点，得到待分类地面点；
[0011]
步骤三，计算每一个待分类地面点的相对高差δh，应用二元高斯混合模型(gmm)进行地物与地面的分类，使用期望最大化算法(em)得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波。
[0012]
所述步骤一的具体过程为：
[0013]
(1)输入机载lidar点云，并根据研究区域点云分布，由用户给定离群点去噪算法(sor) 参数，sor参数包括近邻点个数k以及标准差乘系数m；
[0014]
(2)根据大数定律，遍历每个点云计算到其k近邻点的平均距离，再基于大于m倍标准差为粗差的原则，剔除不满足的点云，得到待处理点云；其中，为第i个点云到其k近邻点的平均距离，knn平均距离标准差标准差n为输入的机载lidar点云数。
[0015]
所述步骤二的具体过程为：
[0016]
(1)先对步骤一得到的待处理点云，由局部最小值法获取初始地面种子点集g
seed
；给定格网大小gw，构建网格索引，根据网格行列号(m,n)计算出网格单元中心的真实平面坐标值(x, y)，再由最小化薄板样条能量函数e
f
的优化原则，选择待插值格网点的n个近邻地面种子点作为控制点；
[0017]
(2)令控制点权重满足插值函数f(x,y)存在二阶偏导数约束条件，根据拉格朗日乘数法，解算出插值函数f(x,y)中的未知系数，进而插值出网格单元中心高程值；
[0018]
(3)给定曲率阈值参数t，由gw和t构建曲面格网金字塔r1，并使用一个3
×
3的均值移动窗口经过每个待处理点云所在格网，得到一个平滑后的曲面格网r
new
，根据r
new
高程值 z
new
(m,n)进而计算平均曲率c(m,n)＝z
new
(m,n)+t；遍历待处理点云，若第j个待处理点云 p
j
(x
j
,y
j
,z
j
)的高程值z
j
大于等于c(m,n)，则将其加入非地面点，否则将其作为第一待分类地面点进入(4)；
[0019]
(4)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔r2，并根据(3)中的方法遍历第一待分类地面点筛选出第二待分类地面点，进入(5)；
[0020]
(5)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔r3，并根据(3)中的方法遍历第二待分类地面点筛选出第三待分类地面点，第一至第三待分类地面点统称待分类地面点。
[0021]
其中，(1)中根据下式换算真实平面坐标值(x，y)：
[0022][0023]
式中，(x
min
，y
min
)为输入的机载lidar点云的最小边界盒的x、y方向上最小坐标值。
[0024]
其中，(2)中插值函数：
[0025][0026]
式中，a1、a2、a3为反映插值曲面总体趋势的系数；w
s
表示与第s个控制点c
s
相关的权重系数；为与薄板样条函数(tps)相关的径向基核函数，(x
s
，y
s
)为c
s
的x、y方向坐标值，n为控制点个数；
[0027]
f(x,y)的二阶偏导数约束为：
[0028][0029]
式中，z
s
为c
s
在z方向的坐标值；
[0030]
根据拉格朗日乘数法，利用以下线性方程组求解f(x,y)中的未知系数：
[0031][0032]
式中，矩阵u的第i行第j列元素(x
i
，y
i
)和(x
j
， y
j
)分别为第i控制点c
i
和第j个控制点c
j
的x、y方向坐标值；的x、y方向坐标值；
[0033]
最后，待插值单元高程值z的计算公式为：z＝f(x,y)。
[0034]
所述步骤三的具体过程为：
[0035]
(1)遍历所有的待分类点，计算每一个待分类点的相对高程，将δh代入二元高斯混合模型(gmm)中，计算高斯混合模型对每个待分类点的相对高程的响应度，应用期望最大
化算法(em)迭代更新计算各分模型参数值α
l
,μ
l
,σ
l2
,l＝1,2；其中，第k个待分类点q
k
的相对高程δh
k
＝h
–
hgrid，h为待分类点的高程值，hgrid为待分类点所在r3的曲面格网单元的高程，α
l
为模型的混合系数且大于等于0，μ
l
是q
k
的相对高程的平均值；σ
l
是q
k
的相对高程的标准差；
[0036]
(2)由于假设地面点、非地面点的相对高程均符合高斯分布，所以根据计算出的各分模型高斯参数值，得到包含隐变量的条件概率估计值，进而使用贝叶斯公式及全概率公式求解用于区分地面点的后验概率p(g|q
k
)，若后验概率p(g|q
k
)大于等于0.5，则认为第k个待分类点q
k
属于地面点，并以之参与数字高程模型dem构建；若后验概率p(g|q
k
)小于0.5，则认为q
k
属于非地面点，并舍弃该点；
[0037]
其中，根据下式更新计算各分模型参数值α
l
,μ
l
,σ
l2
,l＝1,2：
[0038][0039][0040][0041]
式中，δh
k
为待分类地面点q
k
的相对高程；p(g|q
k
)为待分类点q
k
属于地面点的后验概率； p(ng|q
k
)为待分类点q
k
属于非地面点的后验概率；n
c
是待处理点云总数量。
[0042]
其中，(2)中根据下式计算待分类点q
k
属于地面点的后验概率p(g|q
k
)及属于非地面点的后验概率p(ng|q
k
)：
[0043][0044][0045]
式中，p(g)为待分类点q
k
属于地面点的先验概率；p(ng)为待分类点q
k
属于非地面点的先验概率；p(q
k
)由全概率公式计算，p(q
k
)＝p(q
k
|g)p(g)+p(q
k
|ng)p(ng)。
[0046]
每个待分类点q
k
的条件概率：
[0047][0048]
式中，α
l
为模型的混合系数且大于等于0；δh
k
是待处理点q
k
的相对高程；φ()是高
斯密度函数，θ
l
＝(μ
l
,σ
l
)是其参数，μ
l
是待分类点q
k
相对高程的平均值；σ
l
是待分类点q
k
相对高程的标准差。
[0049]
应用上式计算地面点的条件概率p(q
k
|g)，仅用到一组参数(α1,μ1,σ1)：
[0050]
p(q
k
|g)＝α1φ(δh
k
|μ1,σ1)
[0051]
同样地，计算非地面点的条件概率p(q
k
|ng)时用另一组参数(α2,μ2,σ2)：
[0052]
p(q
k
|ng)＝α2φ(δh
k
|μ2,σ2)
[0053]
其中高斯密度函数为：
[0054][0055]
本发明的基于多级曲率约束的高斯混合模型点云滤波方法是一种自适应的滤波方法，其结合了统计方法与曲面插值方法的优点，算法给出了较为通用的一组参数，能够方便有滤波工作需要的用户使用。本发明的滤波方法采用了曲面插值策略，通过薄板样条函数(tps)构建曲面金字塔，避免了因局部地形坡度变化而导致的地形点误分类，通过曲率约束保证了地物点能够被有效区分并被剔除；此外，对每一个测量点进行相对高差的计算，能够有效地克服陡峭地形如悬崖上测量点绝对高程高于平坦地形上的树木测量点绝对高程而导致的误分类问题，进而使用二元高斯混合模型(gmm)通过概率判别的方式完成滤波，模型的求解由期望最大化算法(em)自动地寻优、避免了过多主观判断、及减少了大量的人工调参时间。
附图说明
[0056]
图1是本发明实施例的方法框架图；
[0057]
图2是本发明实施例的统计局外点移除示意图，其中，(a)是移除前示意图，(b)是确定局外点集，(c)是移除后示意图；
[0058]
图3是本发明实施例的多级曲率约束的曲面格网示意图，其中，(a)是选用isprs第三委员会提供的测试样本点云数据所构建的数字表面模型(dsm)，(b)是本发明实施例构建的曲面格网金字塔中第一级曲面格网r1，(c)是本发明实施例构建的曲面格网金字塔中第二级曲面格网r2，(d)是本发明实施例构建的曲面格网金字塔中第三级曲面格网r3；
[0059]
图4是本发明实施例的一个实测点云滤波结果图，其中，(a)是使用原始lidar点云构建的数字表面模型(dsm)，(b)是使用本发明的基于多级曲率约束的高斯混合模型点云滤波方法得到的数字高程模型(dem)；(c)是使用经典数学形态学滤波方法得到的数字高程模型(dem)。
具体实施方式
[0060]
以下结合具体实施例，并参照附图，对本发明作进一步详细说明。
[0061]
本发明根据大数定律，首先计算每个测量点到其k近邻点的平均距离，然后基于正常测量值不会大于m倍标准差的原则，剔除点云中存在的各种异常测量值；其次，对由局部最小值获得的地面种子点，利用薄板样条函数(tps)进行插值，然后剔除与插值曲面高差大于曲率阈值的点。该过程是通过不断减少的插值格网单元与逐渐增加曲率阈值来迭代逼近一个参考曲面；最后，计算每一个测量值的相对高差，应用二元高斯混合模型(gmm)进行地
物与地面的分类，再利用期望最大化算法(em)得到模型参数的最优解，并根据计算得到后验概率判断该点是否属于地物完成点云滤波，以有效地减少人工调参时间，提升机载lidar点云滤波精度与效率。
[0062]
如附图1所示，本发明的基于多级曲率约束的高斯混合模型点云滤波方法包括三部分：
①
根据大数定律确定点云中的待处理点；
②
基于薄板样条插值的多级曲率约束；
③
相对高程计算与高斯混合模型分类。具体实施步骤为：
[0063]
第一步：根据大数定律确定点云中的待处理点。
[0064]
机载lidar点云数据中存在两类噪声点：高位噪声点(high outliers)、低位噪声点(low outliers)，利用本发明中的knn邻域搜索可以检测并剔除这两类噪声。
[0065]
根据大数定律确定点云中的待处理点的示意图如附图2中的(a)至(c)所示，其具体流程如下：
[0066]
(1)输入机载lidar点云，并根据研究区域点云分布，由用户给定离群点去噪算法(sor) 参数，sor参数包括近邻点个数k以及标准差乘系数m；
[0067]
(2)根据大数定律，遍历每个点云计算到其k近邻点的平均距离，再基于大于m倍标准差为粗差的原则，剔除不满足的点云，得到待处理点云；其中，为第i个点云到其k近邻点的平均距离，knn平均距离标准差标准差n为输入的机载lidar点云数。
[0068]
假设待计算点云q
i
(x
i
，y
i
，z
i
)到其最近的k个邻居点p
m
(x
m
，y
m
，z
m
)，(m＝1,2,...,k)的距离为d
im
，则其计算公式如下：
[0069][0070]
则其knn平均距离的计算可以由下式给出：
[0071][0072]
式中，k的值一般由用户指定。
[0073]
再根据可以计算knn平均距离d
mean
及标准差stddev，可以由统计学计算公式给出：
[0074][0075][0076]
式中，n为输入机载lidar点云的总数量。
[0077]
因此，令标准差乘系数m＝3，则待处理点云q
i
的筛选判断条件如下：
[0078]
[0079][0080]
若点满足上式条件约束，则将其加入待处理点集q中；若点不满足上式条件，则将其加入局外点集o并舍弃之。
[0081]
第二步：基于薄板样条插值的多级曲率约束。
[0082]
基于薄板样条插值的多级曲率约束示意图如附图3中的(a)至(d)所示，其过程为：
[0083]
(1)先对步骤一得到的待处理点云，由局部最小值法获取初始地面种子点集g
seed
，给定格网大小gw，构建网格索引，根据网格行列号(m,n)计算出网格单元中心的真实平面坐标值(x, y)，再由最小化薄板样条能量函数e
f
的优化原则，选择待插值格网点的n个近邻地面种子点作为控制点；
[0084]
假设待插值格网的行列号为(m,n)，根据下式换算成相应的真实平面坐标值(x，y)：
[0085][0086]
式中，gw为格网大小；(x
min
，y
min
)为原始点云最小边界盒的x,y方向上最小坐标值。
[0087]
给定点三维空间中的n＝12个控制点c
s
(x
s
，y
s
，z
s
)，薄板样条函数(tps)通过构造函数f(x, y)穿过所有的控制点c
s
并最小化弯曲能量函数e
f
来完成插值。在二维空间中，e
f
被定义为二阶偏导数平方的积分，由下式给出：
[0088][0089]
式中，f为插值函数f(x,y)；(dx,dy)为f(x,y)在x,y方向的微分。
[0090]
(2)令控制点权重满足插值函数f(x,y)存在二阶偏导数约束条件，根据拉格朗日乘数法，解算出插值函数f(x,y)中的未知系数，进而插值出网格单元中心高程值；
[0091]
以待插值单元g
idx
为例，由公式(6)导出其真实平面坐标值(x，y)，其高程值z按下面的方式计算：
[0092]
首先，需要计算f(x,y)中的未知参数，其计算公式由下式给出：
[0093][0094]
式中，a1,a2,a3为反映插值曲面总体趋势的系数；w
s
表示与控制点c
s
相关的权重系数；为与薄板样条函数(tps)相关的径向基核函数，(x
s
，y
s
)为控制点c
s
的x，y方向坐标值；同时，为了满足f(x,y)存在二阶偏导数约束，需要满足下式条件。
[0095][0096]
然后，根据拉格朗日乘数法，利用线性方程组求解f(x,y)中的未知系数，其计算公式如下：
[0097]
[0098]
式中，u为由定义的矩阵，n为控制点个数， (x
i
，y
i
),(x
j
，y
j
)为控制点的x，y方向坐标值；而w为由w
s
构成的列向量；a为由a1,a2,a3构成的列向量；相应地，p和v可以由下式表出：
[0099][0100]
其中，矩阵中的元素来源于一系列控制点集c
s
(x
s
，y
s
，z
s
)，(s＝1,2,
…
,n)，n为控制点个数。
[0101]
最后，待插值单元g
idx
高程值z的计算公式为：
[0102]
z＝f(x,y)
ꢀꢀꢀ
(12)
[0103]
(3)给定曲率阈值参数t，由gw和t构建曲面格网金字塔r1，并使用一个3
×
3的均值移动窗口经过每个待处理点云所在格网，得到一个平滑后的曲面格网r
new
，根据r
new
高程值 z
new
(m,n)进而计算平均曲率c(m,n)＝z
new
(m,n)+t；遍历待处理点云，若点p
j
(x
j
,y
j
,z
j
)的高程值z
j
大于等于c(m,n)，则将其加入非地面点，否则将其作为第一待分类地面点进入(4)。
[0104]
由局部最小值法获取初始地面种子点集g
seed
，然后给定格网大小gw与曲率阈值t，由公式z＝f(x,y)计算出插值曲面格网r1，然后遍历待处理点q
j
，并由公式(6)计算其所在的格网单元行列号(m,n)，采用图像卷积操作，利用一个3
×
3的均值移动窗口处理该格网单元，计算得到一个平滑曲面r
new
，根据r
new
高程值z
new
(m,n)，并按下式计算平均曲率c(m,n)：
[0105]
c(m,n)＝z
new
(m,n)+t
ꢀꢀꢀ
(13)
[0106]
式中，z
new
(m,n)表示行列号(m,n)所在的平滑曲面r
new
格网单元；t为曲率阈值参数。
[0107]
判断待处理点q
j
(x
j
,y
j
,z
j
)的高程值z
j
是否满足条件z
j
<c(m,n)，若是则将其加入待分类地面点集g；否则，将其加入非地面点集ng。
[0108]
(4)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔r2，并根据(3)中的方法遍历第一待分类地面点筛选出第二待分类地面点，进入(5)；
[0109]
(5)令gw＝gw/2且t＝t+0.1，构建曲面格网金字塔r3，并根据(3)中的方法遍历第二待分类地面点筛选出第三待分类地面点，第一至第三待分类地面点统称待分类地面点。
[0110]
构建曲面格网金字塔r1、r2、r3的具体步骤如下：
[0111]
首先，以初始的gw以及t并利用(3)中过程计算得到待分类地面点集g；对其利用局部最小值方法计算更新后的地面种子点集(该步骤目的为了防止线性方程组计算病态问题)，根据更新后的点集建立新的曲面格网，按照公式(13)的曲率约束条件得到新的待分类地面点集g
new
，该过程重复3次结束该级曲面格网r1构建；然后，令gw＝gw/2且t＝t+0.1，构建下一层级的曲面格网r2，对上一步得到的待分类地面点集再次利用局部最小值方法计算更新后的地面种子点集g
new
，根据更新后的点集建立新的曲面格网，按照公式(13)的曲率约束条件得到新的待分类地面点集，该过程重复3次结束该级曲面格网r2构建；最后，令gw＝gw/2 且t＝t+0.1，构建下一层级的曲面格网r3，对上一步得到的待分类地面点集再次利用局部最小值方法计算更新后的地面种子点集g
new
，根据更新后的点集建立新的曲面格网，按照公式 (13)的曲率约束条件得到新的待分类地面点集，该过程重复3次结束该级曲面格网
r3构建。
[0112]
第三步：相对高程计算与高斯混合模型分类。
[0113]
(1)遍历所有的待分类点，计算每一个待分类点的相对高程，将δh代入二元高斯混合模型(gmm)中，计算高斯混合模型对每个待分类点的相对高程的响应度，应用期望最大化算法(em)迭代更新计算各分模型参数值α
l
,μ
l
,σ
l2
,l＝1,2；其中，第k个待分类点q
k
的相对高程δh
k
＝h
–
hgrid，h为待分类点的高程值，hgrid为待分类点所在r3的曲面格网单元的高程；
[0114]
相对高程的计算是为了更显著地区分地面点与非地面点。因为，一般的算法都假设如果一个点高差大于其邻居点的高差，则将其分为非地面点。但是，这一假设当遇到复杂的地形会失效如陡崖等。通过计算每个待处理点q
k
的相对高程δh
k
，可以有效解决这一问题，其公式如下：
[0115]
δh
k
＝h
–
hgrid
ꢀꢀꢀ
(14)
[0116]
式中，h是待处理点q
k
在z方向的坐标值，hgrid是是待处理点q
k
所对应的r3(m,n)的曲面高程值。
[0117]
将δh
k
代入二元高斯混合模型(gmm)中，可以表示每个待处理点q
k
的条件概率：
[0118][0119]
式中，α
l
为模型的混合系数且大于等于0；δh
k
是待处理点q
k
的相对高程；φ()是高斯密度函数，θ
l
＝(μ
l
,σ
l
)是其参数。μ
l
是待处理点q
k
相对高程的平均值；σ
l
是待处理点q
k
相对高程的标准差；其中高斯密度函数为：
[0120][0121]
应用式(15)计算地面点的条件概率p(q
k
|g)，仅用到一组参数(α1,μ1,σ1)：
[0122]
p(q
k
|g)＝α1φ(δh
k
|μ1,σ1)
ꢀꢀꢀ
(17)
[0123]
应用式(15)计算非地面点的条件概率p(q
k
|ng)时用另一组参数(α2,μ2,σ2)：
[0124]
p(q
k
|ng)＝α2φ(δh
k
|μ2,σ2)
ꢀꢀꢀ
(18)
[0125]
由于高斯混合模型中既存在观测变量δh
k
，也存在隐变量g和ng，所以常规的最大似然方法无法被直接用来估计参数。此时，期望最大化算法(em)算法提供了一个求解隐变量近似解的方案。通常，期望最大化算法包含四个主要步骤：
[0126]
i.初始化输入参数，包括α
l
,μ
l
,σ
l
,l＝1,2.
[0127]
首先，我们根据多级曲率约束得到的地面种子点数量以及总待处理点数量n
c
，初始化α1和α2。需要注意的是，α1+α2＝1。而μ1,μ2,σ1,σ2则分别由地面种子点与非地面种子点的相对高程计算。
[0128][0129]
式中，δh
k
的含义与公式(15)一致；n
c
是待处理点云总数量。
[0130]
ii.e步：计算高斯混合模型对每个点的隶属度可以分别表示为p(g|q
k
)和p(ng|q
k
)。
[0131]
iii.m步：计算新一轮迭代中高斯混合模型中的参数(α
l
,μ
l
,σ
l2
,l＝1,2)。
[0132]
迭代中的m步是为了寻找使得e步中隶属度期望最大化的参数值，α
l
,μ
l
,σ
l2
的计算公式如下：
[0133][0134][0135][0136]
式中，δh
k
的含义与公式(15)一致。
[0137]
iv.判断期望最大化算法的收敛条件。
[0138]
期望最大化算法(em)重复(i)-(iv)步骤，直到模型参数(α
l
,μ
l
,σ
l2
,l＝1,2)不再变化为止。
[0139]
(2)由于假设地面点、非地面点的相对高程均符合高斯分布，所以根据计算出的各分模型高斯参数值，得到包含隐变量的条件概率估计值，进而使用贝叶斯公式及全概率公式求解用于区分地面点的后验概率p(g|q
k
)，若后验概率p(g|q
k
)大于等于0.5，则认为第k个待分类点q
k
属于地面点，并以之参与数字高程模型dem构建；若后验概率p(g|q
k
)小于0.5，则认为q
k
属于非地面点，并舍弃该点；
[0140]
由期望最大化算法(em)计算得到模型参数后，首先，直接利用高斯混合模型公式(15-18) 计算地面点的条件概率的估计p(q
k
|g)，及非地面点的条件概率估计p(q
k
|ng)；然后，根据贝叶斯公式，待处理点q
k
的属于地面点的后验概率p(g|q
k
)及属于非地面点的后验概率p(ng |q
k
)的计算公式如下：
[0141][0142][0143]
式中，p(g)代表处理点q
k
属于地面点的先验概率；p(ng)代表处理点q
k
属于非地面点的先验概率；p(q
k
)由全概率公式计算，p(q
k
)＝p(q
k
|g)p(g)+p(q
k
|ng)p(ng)。
[0144]
根据条件p(g|q
k
)≥0.5，若满足则将待分类点q
k
为标记为地面点；否则，则将其分为非地面点，并舍弃之。
[0145]
因此，本发明的基于多级曲率约束的高斯混合模型点云滤波方法是一种基于概率论与数理统计的分类方法，在后验概率的指导下，完成机载lidar点云的滤波。附图4中的(a) 至(c)展示了一个isprs数据集实测点云样本的过滤结果，其中，(a)为是使用原始lidar 点云构建的数字表面模型(dsm)，(b)是使用本发明的基于多级曲率约束的高斯混合模型点云滤波方法得到的数字高程模型(dem)；(c)是使用经典数学形态学滤波方法得到的数字高程模型(dem)。
[0146]
从结果可以看出，经典数学形态学滤波方法得到的数字高程模型(dem)在地形坡度变化剧烈且混合有低矮植被的区域滤波效果不佳，导致许多遗漏误差，一些地面点测量被误分类为非地面点导致错误的dem结果。而本发明的基于多级曲率约束的高斯混合模型点云滤波方法不仅给出了一组较通用的算法参数值，且滤波精度有着明显的改善。
[0147]
以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。