一种基于改进的EM算法的激光雷达波形数据分解的方法

摘要

本发明公开了一种基于改进EM算法的激光雷达波形数据分解的方法，属于机载激光雷达技术领域。通过改进的EM算法即将波形振幅当作一个权值加入到EM算法原始公式的分子和分母上对激光雷达波形进行分析，可以更加详细地了解物体的纵向结构，如表面倾斜、粗糙度、反射率。本发明在后处理过程中让使用者自己提取三维坐标获得高精度的点云结果。本发明使用改进的EM算法得到回波脉冲的位置和宽度，是一种性能可靠、精度较高的波形分解算法。

说明书

技术领域

本发明涉及一种对激光雷达波形数据处理的方法，特别是一种通过改进EM算法实现激光雷达波形数据分解的方法，属于机载激光雷达技术领域。

背景技术

机载激光扫描是一个迅速成长的地形测绘技术。机载激光扫描大多数都重复对地球表面发射红外脉冲，通过光学接收器检测返回传感器的一部分能量。由一个计时器来测量脉冲在激光传感器与地球表面之间的往返时间，通过这个往返时间可以求得传感器到地表的距离，加上经过集成的POS系统和大约2000万～4000万像素的数码相机，使得该技术在一些使用传统摄影测量方法进行地形测绘比较困难的地区，如植被覆盖区、海岸带、岛礁地区等等，发挥着不可替代的作用。

第一个商业机载激光扫描仪只能记录一个后向散射脉冲的时间。如果在激光光斑里只有一个目标那仅仅一个脉冲记录是足够的，在这种情况下返回脉冲可以简单地说明地物。但是对小的激光光斑(0.2-2米)来说，在激光脉冲的传播路径里可能有若干地物产生各自的后向散射脉冲。已经有更加先进的激光扫描仪能够记录多次回波，但这种工作方式使用户无法获得任何与有关设备本身相关的一些信息的，如如何根据回波信号定位，电子设备和地物结构对获取的回波信号的形状和大小有何影响，回波信号如何被量化成几次离散脉冲信号等。另一方面，机载激光雷达的光斑设计得越来越小，这样每束激光产生多次回波的数量也进一步增加，这对脉冲探测和量化方法提出了更大的挑战。曾有文献指出，探测和量化方法不同，可以引起成果的误差。

解决上述问题的途径是将发射信号和回波信号均以很小的采样间隔进行采样并记录，而不仅仅是记录若干次离散的回波信号。这样的采样记录方式即所谓的全波形数字化记录，这种类型的激光雷达系统称为waveform-digitizingLIDAR。事实上早在上世纪90年代，NASA开发的一些机载激光雷达系统已经具备这个能力，比如SLICER(Scanning Lidar Imager of Canopies by Echo Recover)和LVIS(Laser Vegetation Imaging Sensor)。还有一些星载激光雷达系统如GLAS也具有全波形记录能力。但是与硬件发展相比，波形数据的分析和处理方法的研究却相对滞后。这一方面是由于波形数据尚未真正普及，另一方面由于波形数据的分析与具体应用关系密切，分析方法要面向应用。

发明内容

本发明的目的就是提供一种处理激光雷达波形数据的方法，根据激光雷达波形数据符合高斯分布的特征，利用改进的EM(Expectation-Maximization)算法高斯分解激光雷达波形数据，生产出高质量的点云数据并得到森林参数。

实现本发明目的采用的技术方案：一种基于改进EM算法的激光雷达波形数据分解的方法，包括以下步骤：

(1)依据雷达等式的推导，用高斯函数模拟激光雷达的波形数据；

(2)通过使用去噪阈值和平滑算法去除噪音；

(3)在去噪后会出现毛刺现象，将毛刺作为噪音再进行一次去噪处理；

(4)在去除噪音后，通过波形维护来优化波形；

(5)对完全去除噪音后的激光雷达波形数据求梯度算子，把检测到的若干激光雷达波形极大值作为高斯函数的初始值；

(6)通过EM算法对高斯函数参数做最大相似评估，将波形振幅当作一个权值加入到EM算法原始公式的分子和分母上；

(7)用改进的EM算法高斯分解波形数据，由于激光雷达波形数据会产生波形重叠的情况，运用最小距离法通过计算高斯函数期望值之间的差值来决定重叠波形的高斯分解，所述改进的EM算法即为将波形振幅当作一个权值加入到EM算法原始公式的分子和分母上；

(8)提取分解激光雷达波形得到的高斯函数，高斯函数期望值为激光雷达波形的位置，高斯函数的均方差为激光雷达波形的宽度；

(9)通过应用高斯函数期望值求解激光雷达点云数据三维坐标，均方差求解森林参数等地物属性。

在步骤(2)和(3)中，去噪处理包括以下步骤：

(1)通过求得发射脉冲波形的背景噪音平均值和均方差来确定返回脉冲波形的去噪阈值，在返回脉冲波形里将小于去噪阈值的脉冲振幅赋值为零；

(2)应用平滑算法把波形曲线上的噪音抖动部分去掉。

在步骤(4)中，波形维护包括以下步骤：

(1)检测去噪后每一段波形数据的位置；

(2)在经过平滑的原始波形数据上找到相应的波形段，并把这些波形段两端的数据恢复到去噪后的波形数据里。

本发明通过用高斯分解的方式对激光雷达单个光斑范围里的地物反射回波脉冲进行量化描述，实现对地物位置、反射率、粗糙度等一系列属性提取，且利用改进的EM算法高斯分解激光雷达波形数据，生产出高质量的点云数据并得到森林参数。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是原始波形数据。

图2是带有毛刺的去噪波形。

图3是平滑后的去噪波形。

图4是预处理后的波形。

图5是维护前的波形。

图6是维护后的波形。

图7是波形的一阶导数曲线。

图8是设定初值后的波形图。

图9是基于EM算法波形高斯分解的初始流程图。

图10是基于EM算法的波形模拟结果。

图11是改进后的EM算法模拟的波形。

图12是SLICER系统提取的地表位置。

图13是改进后的EM算法模拟的SLICER波形。

图14是基于改进后的EM算法波形高斯分解的最终流程图。

具体实施方式

本发明提供的基于改进的EM算法的激光雷达波形数据处理的方法包括以下步骤：

(1)根据雷达等式的推导，用高斯函数模拟激光雷达的波形数据是一个最佳方案，激光雷达原始波形数据如图1所示，横坐标为波形数字化的采样数，纵坐标为采样的振幅值。

(2)去除波形数据的噪音。波形在采集的过程中由于多方面的原因会存在噪音，所以为一个抖动的曲线，图1中振幅很小且抖动的部位即为噪音。先设定一个阈值来去除噪音，去噪分两步来处理：第一步通过求得发射脉冲波形的背景噪音平均值和均方差来确定返回脉冲波形的去噪阈值σ_noise，在返回脉冲波形里小于去噪阈值σ_noise的脉冲振幅赋值为零。由于波形的噪音振幅随机性非常的大，求取的去噪阈值σ_noise无法把所有的噪音都去掉，造成去噪后的波形有毛刺，如图2所示。第二步应用平滑算法把波形曲线上的噪音抖动部分去掉，如图3所示，去除噪音后的波形曲线是一条流畅平滑的信号曲线，这里可以看到有效波形部分的噪音通过平滑算法被去掉，毛刺部分也因为平滑算法而被压下去，但毛刺依然没有去除。

(3)由于有些波形数据的噪音部分抖动得很厉害，在去噪后会出现很多毛刺现象，检测出毛刺，并将它们作为噪音再进行一次去噪处理。再用去噪阈值σ_noise做一次去噪处理，最终得到的结果如图4所示。

(4)如图5所示，由于有些地物的反射率很低或者由于遮挡打在地物上的能量很弱使其反射很小的能量，在去除噪音的时候使有效波形损失过大不利于后续处理。对去噪后的波形进行维护分两步进行：第一步检测去噪后每一段波形数据的位置，第二步在经过平滑的原始波形数据上找到这些相应的波形段，并把这些波形段两端一定范围的数据恢复到去噪后的波形数据里，图6为维护后的波形。图5和图6中的横线为去噪阈值σ_noise。

(5)使用EM算法之前，μ_j的初始值为波形的局部极大值，如图6所示的位置，因此可运用梯度算子对波形求一阶导数，图7是图6的一阶导数曲线，左边大于0右边小于0的采样点就是局部极大值所在的位置，并由极大值的个数确定初始总的高斯函数数量。p_j的初始值是让所有的j个高斯分支有相同的权，在本实施例中σ_j的初始值设置为1。如图8所示，线1是预处理后的波形数据，线2是设定初始值后的高斯函数曲线。

(6)由于回波波形可以看作是若干高斯函数的叠加，波形的分解即为混合高斯分布，通过EM算法完成高斯混合密度的参数估计。

EM算法计算p_j、μ_j、σ_j最优值的原始公式为：

$Q_{\mathrm{ij}}=\frac{p_j{f}_j\left(x_i\right)}{{\Sigma }_{j=1}^k{p}_j{f}_j\left(x_i\right)}---\left(1\right)$

$p_j=\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}}{n}---\left(2\right)$

${\mu }_j=\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}i}{p_j\times n}---\left(3\right)$

${\sigma }_j=\sqrt{\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}{(i-{\mu }_j)}^2}{p_j\times n}}---\left(4\right)$

其中Q_ij是采样i属于高斯函数分支j的权重，f_j(x)是高斯概率密度函数，p_j是f_j(x)的权，μ_j是高斯函数的期望值，σ_j是高斯函数的标准差，通过求得的μ_j、σ_j初始值得到初始高斯函数乃f_j(x)～N(μ_j，σ_j²)，波形由公式$f\left(x\right)=\underset{j=1}{\overset{k}{\Sigma }}{p}_j\times f_i\left(x\right)$来描述。通过迭代公式(1)～(4)不断地调节p_j、μ_j、σ_j这三个参数得到最优值。波形数据分解的初始流程图如图9所示。得到波形分解的结果图如图10所示，左边的高斯模拟波形是合理的，但右边存在很大的偏差。

波形由一系列简单的高斯分布组成，则混合分布的数学表示式(即采样后得到的波形)可以由下式表示：

$f\left(x\right)=\underset{j=1}{\overset{k}{\Sigma }}{p}_j\times f_i\left(x\right),$$f_i\left(x\right)\in N({\mu }_j,{\sigma }_j^2)$

其中k是高斯函数的数量，f_j(x)是高斯概率密度函数，p_j是f_j(x)的权，表示该分布在混合分布中占的比重，满足：$0<p_j<1,\underset{j=1}{\overset{k}{\Sigma }}{p}_j=1,$μ_j是斯函数期望值，σ_j是高斯函数的标准差。对每个部分j，计算得到的μ_j表示回波的位置，σ_j表示回波的宽度。参数p_j、μ_j、σ_j均可以由EM算法估计得到。先对数据进行预处理，并对公式(1)～(4)的参数定义初始值。由于初始值是预测的，f(x)也是依据初始值μ_j、σ_j得到的高斯函数，然后用f(x)来计算Q_ij，但在后来的公式(1)-(4)中还是仅仅用预测的Q_ij来计算，而不考虑实际的振幅值来拟合波形数据得不到一个理想结果。所以从物理意义上来说在公式(1)～(4)里加上一个振幅值N_i。以下为公式推导：

将公式(2)代入到公式(3)～(4)中得到：

${\mu }_j=\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}i}{n\times p_j}=\frac{n\times {\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}i}{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}\times n}=\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}i}{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}}---\left(5\right)$

${\sigma }_j=\sqrt{\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}{(i-{\mu }_j)}^2}{n\times p_j}}=\sqrt{\frac{n\times {\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}{(i-{\mu }_j)}^2}{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}\times n}}=\sqrt{\frac{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}{(i-{\mu }_j)}^2}{{\Sigma }_{i=1}^n{Q}_{\mathrm{ij}}}}---\left(6\right)$

现在同时在分母和分子加入振幅N_i，N_i相当于一个权值，用来约束p_j、μ_j、σ_j值，得到公式如下：

${\mu }_j=\frac{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}i}{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}}---\left(7\right)$

${\sigma }_j=\sqrt{\frac{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}{(i-{\mu }_j)}^2}{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}}}---\left(8\right)$

把公式(7)-(8)还原为EM算法的形式：

$Q_{\mathrm{ij}}=\frac{p_j{f}_j\left(i\right)}{{\Sigma }_{j=1}^k{p}_j{f}_j\left(i\right)}---\left(9\right)$

$p_j=\frac{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}}{n\times {\Sigma }_{i=1}^n{N}_i}---\left(10\right)$

${\mu }_j=\frac{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}i}{n\times p_j\times {\Sigma }_{i=1}^n{N}_i}---\left(11\right)$

${\sigma }_j=\sqrt{\frac{{\Sigma }_{i=1}^n{N}_i{Q}_{\mathrm{ij}}{(1-{\mu }_j)}^2}{n\times p_j\times {\Sigma }_{i=1}^n{N}_i}}---\left(12\right)$

n是波形中采样的数量，N_i是第i次采样的振幅。图11为用改进的EM算法得到的结果，可以看出其很好地模拟了波形数据。

将改进的EM算法得到的波形参数与SLICER系统得到的波形参数进行比较：

SLICER提供了一个峰值提取算法，如图12所示，三条竖线为系统提取的Grstart、Grpeak、GrEnd三个峰值，分别代表地表返回的开始、峰值和末端。在Grpeak、GrEnd之间有一个波峰，这个波峰是真正的地表返回，这束激光产生了两次明显的回波：首先其打在了一个浓密的低矮植被上产生一个振幅较大的回波，然后打在地表。由于大部分激光能量被植被反射，所以地表反射的回波较小。又由于植被与地表距离较近，使两次回波发生了部分重叠。由此可见系统提取的地表返回位置不够精确。

用改进的EM算法得到的模拟波形如图13所示，虚线为EM算法模拟的两次回波，两条纵向直线标明了两次回波波峰的位置，μ₁＝57、σ₁＝7.3、μ₂＝83，σ₂＝7.3。将图13与图12进行对比，图13得到的结果精度更高，正确性更好。所以利用改进的EM算法波形分解可以获得更高的精度和更好的正确性。

(7)在找到最优值后，采用min dist(k)＝min(μ_j+1-μ_j)方法，确定得到的k值是否是最佳的，即用k个高斯函数描述波形是否合适。图14为是基于改进后的EM算法波形高斯分解的最终流程图。

(8)得到波形高斯分解的期望值后，依据波形数据起始点的坐标可以计算点云数据，波形分解得到的点云数据精度要比系统RIEGL LiteMapper5600的要高，波形分解的点云数据也更加有层次感，奇异点的数量也明显减少，其可以为DTM生产、城市建模等提供优质的数据源。高斯分解得到的点云无论是数量还是质量都优于系统RIEGLLiteMapper5600提取的点云数据。

(9)在森林覆盖地区，用高斯分解的最后一次回波到反射波形开始处之间的距离来表示树高；用激光打在树冠上的脉冲总量来描述冠层蓄积量；单次回波用高斯分解检测，将检测出的波形数据用公式计算地表反射率；将检测出多次回波的数据用公式计算植被反射率，由于影响反射率的因数很多，用地表平均反射率来计算植被平均反射率会引起误差，但这个平均反射率在局部地区还是相当有代表性的，其结果可以用于分类；应用公式来计算森林郁闭度。