CN104635223A - 基于经验模式分解与分数阶傅里叶变换的激光回波去噪方法 - Google Patents

Abstract

本发明提供了一种基于经验模式分解与分数阶傅里叶变换的激光雷达回波去噪算法，包括经验模式分解、选取高频及过渡固有模态函数、分数阶傅里叶变换、实际信号分数阶傅里叶变换频谱分析、带通滤波、逆分数阶傅里叶变换以及信号重构。通过经验模式分解将信号分解为不同频带尺度的时域信号分量，选取特定的固有模态函数(Intrinsic Mode Function，IMF)，对其作分数阶傅里叶变换，滤波后重构信号，提高去噪效果。

Description

基于经验模式分解与分数阶傅里叶变换的激光回波去噪方法

技术领域

本发明属于激光雷达数据处理技术领域，具体是指一种结合了经验模式分解、频谱分析以及分数阶傅里叶变换的激光雷达回波去噪方法。适用于激光雷达回波数据的预处理。

背景技术

激光雷达回波信号中不可避免的会存在各种噪声，噪声的存在影响激光数据的处理和解释。噪声的压制与消除有着多种多样的方法，其中典型的代表就是小波域方法，小波滤波的机理是基于信号与噪声的小波系数在尺度上的不同性质，采用相应准则对含噪小波系数进行处理，从而达到去噪目的，是近二十几年来时频分析研究的热点。该方法的基本步骤如下：首先选择合适的小波基、确定分解层次、对含有噪声的原始信号作小波变换，然后选取阈值函数，利用门限处理相应的小波系数，获得新的小波系数，最后进行小波逆变换，获得去噪后的信号。然而小波域变换在本质上仍然是一种线性变换，并且对于小波基的选择具有依赖性，不具有自适应性。而激光雷达信号是一种典型的非线性非平稳信号。

Huang等人提出了一种新的时域信号处理方法：经验模式分解(Empirical Mode Decomposition,EMD)，EMD的最大优点是可以自适应的从数据中得到基函数，克服了小波变换中要选取合适的小波基的困难，根据信号时间尺度的不同，EMD将信号分解为几个不同频带尺度的时域信号分量，从瞬态尺度(高频模态)到粗糙尺度(低频模态)。实际应用中，噪声往往处于高频部分，然而，简单的舍弃高频分量会导致有效信号的损失，而舍弃过少的高频分量会导致信号去噪效果不明显。

发明内容

本发明的目的是基于激光雷达系统，开发一种具有良好的去噪特性的基于经验模式分解及分数阶傅里叶变换的激光雷达回波去噪方法。

本发明的技术方案是：一种基于经验模式分解及分数阶傅里叶变换的激光雷达回波去噪方法，包括以下步骤：

(1)依据实际所用激光系统，获取激光雷达回波数据信号；

(2)对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数，并确定其中的高频及过渡固有模态函数；

(3)对系统实际打出的不含噪声的理想激光回波信号作分数阶傅里叶变换，分析这一阶上的信号频谱；

(4)对步骤(2)中所得的高频及过渡固有模态函数作分数阶傅里叶变换，并根据所得理想激光回波信号的频谱分析，对进行分数阶傅里叶变换后的高频及过渡固有模态函数进行带通滤波；

(5)对带通滤波后的信号进行逆分数阶傅里叶变换；

(6)将步骤(5)中处理后的固有模态函数分量与步骤(2)中除高频及过渡固有模态函数外的剩余未处理的固有模态函数分量重构，得到去噪后的信号。

进一步的，步骤(2)中对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数，若分解得到的第一固有模态函数中噪声含量接近全部时，则直接删去，再在剩余的固有模态函数中确定高频及过渡固有模态函数。

进一步的，步骤(2)对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数；具体如下：

对于一个激光雷达回波数据信号x(t)，假设原始残差信号为r₀(t)＝x(t)，通过EMD获得固有模态函数的步骤为：

②识别出信号x(t)的所有极大值点和极小值点，分别拟合出信号x(t)的上、下包络线e_upp(t)、e_low(t)，计算上下包络线的平均值m₁(t)：

$m_1\left(t\right)=\frac{e_{\mathrm{upp}}\left(t\right)+e_{\mathrm{low}}\left(t\right)}{2}$

②将x(t)减去m₁(t)得到h₁(t)，将h₁(t)视为新的x(t)，重复步骤①，激光雷达回波数据信号x(t)经过k次筛选，直到h_k(t)是基本IMF分量，即从原始数据中处理得到第一个固有模态函数分量h_1,k(t)，它包含原始信号中最短的周期分量；

③定义c₁(t)，c₁(t)＝h_1,k(t)，从原始信号中分离出c₁(t)，得到：

r₁(t)＝x(t)-c₁(t)；

④因剩余部分r₁(t)仍然包含较长周期分量的信息，所以r₁(t)仍被当作新的x(t)按以上相同的处理过程来处理；该处理过程可对所有的接下来的剩余分量r_j(t)进行处理，得到如下结果：

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ r_2\left(t\right)-c_3\left(t\right)=r_3\left(t\right)\\ .\\ .\\ .\\ r_{N-1}\left(t\right)-c_N\left(t\right)=r_N\left(t\right)\end{array}\right\}$

⑤x(t)最终表示为：

$x\left(t\right)=\underset{1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_N\left(t\right),i=1\·\·\·N.$

本发明的有益效果是：将EMD与分数阶傅里叶变换(Fractional Fourier Transform,FRFT)相结合，对数据进行处理，结果表明该算法具有良好的去噪特性；与传统的经验模式分解舍弃第一个高频分量对比，信噪比的提升效果很大(信噪比由15.8033dB提升20.0741dB，而对于信噪比更低的信号，信噪比的提升效果更大)。

附图说明

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

图1是本发明的技术流程图。

图2是实际激光雷达回波信号。

图3是对信号作经验模式分解后得到的前八层IMF分量，自上而下分别为IMF1到IMF8。

图4是对理想信号作0.5阶分数阶傅里叶变换后的结果，观察其频谱特征。

图5(a)是选取IMF2分量分量后对其作0.5阶分数阶傅里叶变换后的结果，并对其作逆分数阶傅里叶变换后的结果。

图5(b)是选取IMF3分量分量后对其作0.5阶分数阶傅里叶变换后的结果，并对其作逆分数阶傅里叶变换后的结果。

图6是对处理后的IMF分量与剩余IMF分量进行重构后的去噪结果。

具体实施方式

下面结合附图对本发明的技术实施方案进行详细说明。

如图1所示，一种基于经验模式分解与分数阶傅里叶变换的激光回波去噪方法，包括以下步骤：

(1)依据实际所用激光系统，获取激光雷达回波数据信号如图2所示，横坐标为波形数据的采样间隔(单位:ns)，纵坐标为幅度值。

(2)对回波数据信号x(t)作经验模式分解，得到若干个不同层次的固有模态函数，确定高频及过渡分量。

经验模式分解的假设条件是所要分解的信号是由许多不同的固有模态函数(IMF)叠加而成的，每个IMF代表一个振动模式，各IMF分量包含了原信号的不同时间尺度的局部特征信号，IMF要满足以下两个条件：

a、整个IMF中零点数与极点数相等或至多相差1；

b、IMF上任意一点由局部极大值点确定的上包络线和由局部极小值点确定的下包络线的均值都为零，即信号关于时间轴局部对称。

对于一个激光雷达回波数据信号x(t)，假设原始残差信号为r₀(t)＝x(t)，通过EMD获得IMF的步骤为：

③识别出信号x(t)的所有极大值点和极小值点，分别拟合出信号x(t)的上、下包络线e_upp(t)、e_low(t)，计算上下包络线的平均值m₁(t)：

$m_1\left(t\right)=\frac{e_{\mathrm{upp}}\left(t\right)+e_{\mathrm{low}}\left(t\right)}{2}$

②将x(t)减去m₁(t)得到h₁(t)，将h₁(t)视为新的x(t)，重复步骤①，激光雷达回波数据信号x(t)经过k次筛选，直到h_k(t)是基本IMF分量，即从原始数据中处理得到第一个固有模态函数分量h_1,k(t)，它应包含原始信号中最短的周期分量；

③定义c₁(t)，c₁(t)＝h_1,k(t)，从原始信号中分离出c₁(t)，得到：

r₁(t)＝x(t)-c₁(t)

④因剩余部分r₁(t)仍然包含较长周期分量的信息，所以r₁(t)仍被当作新的x(t)按以上相同的处理过程来处理。该处理过程可对所有的接下来的剩余分量r_j(t)进行处理，得到如下结果：

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ r_2\left(t\right)-c_3\left(t\right)=r_3\left(t\right)\\ .\\ .\\ .\\ r_{N-1}\left(t\right)-c_N\left(t\right)=r_N\left(t\right)\end{array}\right\}$

⑤x(t)最终表示为

$x\left(t\right)=\underset{1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_N\left(t\right),i=1\·\·\·N.$

综上所述，EMD的分解过程实际上是一个筛分过程，在筛分过程中，不仅消除了模态波形的叠加，而且使波形轮廓更加对称。EMD方法从特征时间尺度出发，首先将信号中特征时间尺度最小的模态分离出来，然后分离特征时间尺度较大的模态，最后分离特征时间尺度最大的模态。分解结果如图3所示。因IMF1中几乎全部为噪声，所以直接删去。

(3)对系统实际打出的不含噪声的理想激光回波信号作分数阶傅里叶变换。

分数阶傅里叶变换作为一种较新的时频分析工具，是一种广义化的傅里叶变换，它既与经典的傅里叶变换有着联系，又提供了某些其不具有的特点。分数阶傅里叶变换将信号分解到一系列同一调频率不同起始频率的线性调频正交基上，在同一的时频域上进行信号处理，灵活性更强，适合于非平稳信号的处理。

FRFT的基本定义为 $X_p\left(u\right)={\∫}_{-\∞}^{+\∞}{\tilde{K}}_p(u,t)x\left(t\right)\mathrm{dt}.$ 其中，

${\tilde{K}}_p(u,t)=\left\{\begin{array}{c}{A}_{\mathrm{\α}}\exp \left[\mathrm{j\π}(u^2\cot \mathrm{\α}-2\mathrm{ut}\csc \mathrm{\α}+t^2\cot \mathrm{\α})\right],\mathrm{\α}\≠\mathrm{n\π}\\ \mathrm{\δ}(u-t),p=4n(\mathrm{\α}=2\mathrm{n\π})\\ \mathrm{\δ}(u+t),p=4n\±2(\mathrm{\α}=(2n\±1)\mathrm{\π})\end{array}\right.$ 为FRFT的核函数，x(t)为待转换的时域信号，X_p(u)为分数阶傅里叶变换后的频域信号。

α＝pπ/2,p≠2n，p为分数阶次，n为整数。当p＝1时，A_α＝1,α＝π/2,且为普通的傅里叶变换，因而分数阶傅里叶变换可以看做是广义的傅里叶变换变换的一种推广。

本发明中选取分数阶次p＝0.5，理想信号的0.5阶分数阶傅里叶变换频谱如图4所示。

(4)选取高频及过渡IMF分量，本发明中选取两个IMF分量：IMF2以及IMF3分量，对其作0.5阶分数阶傅里叶变换，其结果如图5所示。

(5)根据理想激光回波信号的0.5阶分数阶傅里叶变换频谱特征，对IMF2以及IMF3的分数阶傅里叶变换频谱进行带通滤波，并作逆分数阶傅里叶变换。将处理后的IMF2以及IMF3与剩余的未处理的IMF4-8分量进行重构，得到去噪后的信号，结果如图6所示。与传统的经验模式分解舍弃第一个高频分量对比，信噪比由15.8033dB提升20.0741dB，而对于信噪比更低的信号，信噪比的提升效果将更大。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明。本发明所属技术领域中具有通常知识者，在不脱离本发明的精神和范围内，当可作各种的更动与润饰。因此，本发明的保护范围当视权利要求书所界定者为准。

Claims (3)

1.基于经验模式分解及分数阶傅里叶变换的激光雷达回波去噪方法，其特征在于：包括以下步骤：

(1)依据实际所用激光系统，获取激光雷达回波数据信号；

(2)对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数，并确定其中的高频及过渡固有模态函数；

(3)对系统实际打出的不含噪声的理想激光回波信号作分数阶傅里叶变换，分析这一阶上的信号频谱；

(4)对步骤(2)中所得的高频及过渡固有模态函数作分数阶傅里叶变换，并根据所得理想激光回波信号的频谱分析，对进行分数阶傅里叶变换后的高频及过渡固有模态函数进行带通滤波；

(5)对带通滤波后的信号进行逆分数阶傅里叶变换；

(6)将步骤(5)中处理后的固有模态函数分量与步骤(2)中除高频及过渡固有模态函数外的剩余未处理的固有模态函数分量重构，得到去噪后的信号。

2.根据权利要求1所述的基于经验模式分解及分数阶傅里叶变换的激光雷达回波去噪方法，其特征在于：步骤(2)中对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数，若分解得到的第一固有模态函数中噪声含量接近全部时，则直接删去，再在剩余的固有模态函数中确定高频及过渡固有模态函数。

3.根据权利要求1所述的基于经验模式分解及分数阶傅里叶变换的激光雷达回波去噪方法，其特征在于：步骤(2)对获取的激光雷达回波数据信号作经验模式分解，得到若干个不同层次的固有模态函数；具体如下：

对于一个激光雷达回波数据信号x(t)，假设原始残差信号为r₀(t)＝x(t)，通过EMD获得固有模态函数的步骤为：

①识别出信号x(t)的所有极大值点和极小值点，分别拟合出信号x(t)的上、下包络线e_upp(t)、e_low(t)，计算上下包络线的平均值m₁(t)：

$m_1\left(t\right)=\frac{e_{\mathrm{upp}}\left(t\right)+e_{\mathrm{low}}\left(t\right)}{2}$

②将x(t)减去m₁(t)得到h₁(t)，将h₁(t)视为新的x(t)，重复步骤①，激光雷达回波数据信号x(t)经过k次筛选，直到h_k(t)是基本IMF分量，即从原始数据中处理得到第一个固有模态函数分量h_1,k(t)，它包含原始信号中最短的周期分量；

③定义c₁(t)，c₁(t)＝h_1,k(t)，从原始信号中分离出c₁(t)，得到：

r₁(t)＝x(t)-c₁(t)；

④因剩余部分r₁(t)仍然包含较长周期分量的信息，所以r₁(t)仍被当作新的x(t)按以上相同的处理过程来处理；该处理过程可对所有的接下来的剩余分量r_j(t)进行处理，得到如下结果：

$\left.\begin{array}{c}{r}_1\left(t\right)-c_2\left(t\right)=r_2\left(t\right)\\ r_2\left(t\right)-c_3\left(t\right)=r_3\left(t\right)\\ .\\ .\\ .\\ r_{N-1}\left(t\right)-c_N\left(t\right)=r_N\left(t\right)\end{array}\right\}$

⑤x(t)最终表示为：

$x\left(t\right)=\underset{1}{\overset{N}{\mathrm{\Σ}}}{c}_i\left(t\right)+r_N\left(t\right),i=1...N.$