一种保形LiDAR波形去噪方法及系统

摘要

本发明提供一种保形LiDAR波形去噪方法及系统，包括迭代进行去噪过程，每次去噪依次对每个点，先通过对邻域计算得出该采样点的估算值，将原始值与其进行比较和判定，若差值大于所设定的阈值，则将估算值赋给采样点，否则保持原始值；估算值的获取方式为，首先判断当前操作数为奇数还是偶数，当操作数为奇数时和当操作数为偶数时进行相反的计算；在当前的操作数达到预设的迭代次数时结束流程，得到去噪结果。本发明技术方案能够有效去除LiDAR回波信号中的随机噪声；能够在去噪同时保持波形的形状。应用本发明技术方案对LiDAR回波信号进行去噪可以有效提高信噪比，利于后续处理。

说明书

技术领域

本发明属于数字信号去噪处理领域，特别涉及一种保形LiDAR波形去噪方法及系统。

背景技术

目前，方海涛等在频率域通过采用非线性软阈值技术的小波系数去噪方法去除了信号中 的噪声，并保持信号的特征(Optics Communication 233卷2004年)。除了在频率域，还可以 在空间域对LiDAR回波信号进行去噪。孙冰玉等将LiDAR信号去噪问题视为一个函数回归 问题，利用最小二乘支持向量机对LiDAR信号进行去噪(Process Letters 12卷2期2005年)。 这些技术都很好地实现了对LiDAR回波信号的去噪，但是频率域的去噪若用于LiDAR成千 上万次回波的去噪，消耗的时间空间太大，传统的图像空间域处理滤波方法不需要波形的先 验知识，且方法简单适用于快速处理大量波形数据，但是这些方法在去噪的同时，可能会使 波形发生畸变，如波形的峰值大小、波宽等。

发明内容

针对以上技术问题，本发明从空间域出发，在邻域法的基础上设计了一个改进的邻域光 顺去噪技术方案，优化了去噪效果，减少了波形的变形。

为达到上述目的，本发明采用如下的技术方案：

一种保形LiDAR波形去噪方法，包括如下步骤：

步骤1，对于一维信号x，设包括n个点x₁,…,x_n，标记x_i为一维信号x中第i个点，令操作数 m＝1；

步骤2，进行第m次去噪过程，包括依次对每个点x_i，先通过对邻域计算得出该点的估算 值将x_i原始值与其进行比较和判定，若差值大于所设定的阈值，则将赋给 点x_i，否则保持原始值；i的取值为2,3,…,n-1，x₁和x_n保持不变；

估算值的获取方式为，首先判断当前操作数为奇数还是偶数，

当操作数为奇数时进行以下计算，

x′_i＝x_i+λΔx_i，当i为奇数时

x′_i＝x_i+μΔx_i，当i为偶数时

当操作数为偶数时进行如下计算，

x′_i＝x_i+λΔx_i，当i为偶数时

x′_i＝x_i+μΔx_i，当i为奇数时

其中，变量λ是一个大于0小于1的尺度因子，μ是一个负尺 度因子，且μ<-λ；将计算所得x′_i作为估算值

步骤3，判断当前的操作数m是否达到预设的迭代次数，是则结束流程，得到去噪结果； 否则令m＝m+1，返回步骤2进行迭代，执行下一次去噪过程。

本发明还相应提供一种保形LiDAR波形去噪系统，包括如下模块：

初始化模块，用于对于一维信号x，设包括n个点x₁,…,x_n，标记x_i为一维信号x中第i个点， 令操作数m＝1；

去噪模块，用于进行第m次去噪过程，包括依次对每个点x_i，先通过对邻域计算得出该点 的估算值将x_i原始值与其进行比较和判定，若大于所设定的阈值，则将赋 给点x_i，否则保持原始值；i的取值为2,3,…,n-1，x₁和x_n保持不变；

估算值的获取方式为，首先判断当前操作数为奇数还是偶数，

当操作数为奇数时进行以下计算，

x′_i＝x_i+λΔx_i，当i为奇数时

x′_i＝x_i+μΔx_i，当i为偶数时

当操作数为偶数时进行如下计算，

x′_i＝x_i+λΔx_i，当i为偶数时

x′_i＝x_i+μΔx_i，当i为奇数时

其中，变量λ是一个大于0小于1的尺度因子，μ是一个负尺 度因子，且μ<-λ；将计算所得x′_i作为估算值

迭代判断模块，用于判断当前的操作数m是否达到预设的迭代次数，是则结束工作，得 到去噪结果；否则令m＝m+1，命令去噪模块进行迭代工作，执行下一次去噪过程。

本发明技术方案能够有效去除LiDAR回波信号中的随机噪声；能够在去噪同时保持波形 的形状。应用本发明技术方案对LiDAR回波信号进行去噪可以有效提高信噪比，利于后续处 理。

具体实施方式

传统λ/μ去噪方法的基本原理是对采样点及其邻域进行适当操作，使整个波形平滑，曲 率变化连续。因此一个点是否为噪声点，是可以通过看它是否与其邻域点连续来判断。若为 含噪声点，则应将其去掉，并对其邻域值进行合适的计算，重新给它赋予一个近似真实的值。 通常情况下，去噪过程需迭代进行，本领域技术人员可自行预设合适的操作执行次数。

每次操作，依次对每个点，先通过对邻域计算得出该采样点的估算值将原始值x_i与 其进行比较和判定。若差值大于所设定的阈值(本领域技术人员可自行预设阈值)， 则将赋给点x_i，否则保持原始值。

对于一个一维信号x，标记x_i为一维信号x中第i个点，设包括n个点x₁,…,x_n，即i的取值为 1,2,…,n，可将其表示成一个列向量(x₁,…,x_n)^t，高斯滤波的最简单形式可表述为式(1)：

x′_i＝x_i+λΔx_i  (1)

其中，x′_i为x_i高斯滤波后的结果，i的取值为2,3,…,n-1，x₁和x_n保持不变；变量 λ是一个大于0小于1的尺度因子。该方程可以写成矩阵形式， 如式(2)：

x′＝(I-λK)x  (2)

其中，x′为一维信号x高斯滤波后的结果，矩阵I为单位矩阵，矩阵K如式(3)：

矩阵K为n行n列，n为一维信号的采样点数，矩阵中除了主对角线上的值为2，其余值都 是-1。

高斯滤波会导致信号边缘的收缩，所以高斯核函数不能作为低通滤波的核，于是采用一 个非收缩的关于矩阵K的函数f(K)来代替矩阵I-λK，此时：

x′＝f(K)x  (4)

如果迭代N次，则结果如式(5)：

x^N＝f(K)^N  (5)

其中，x^N为一维信号x经N次迭代进行高斯滤波后的结果。

由于矩阵K是对称的，它具有实特征值和特征向量。设k_i为矩阵K的第i个实特征值，k_i对 应的特征向量为u_i，矩阵K的实特征值为0≤k₁≤k₂≤...≤k_n≤m，m为常数，k₁,...,k_n对应的特 征向量分别为u₁,...,u_n，式(4)可表示为式(6)：

$x^{\prime }=f\left(K\right)x={\Sigma }_{i=1}^n{\xi }_{i}f\left(k_i\right)u_i---\left(6\right)$

其中，ξ_i为矩阵K的矩阵函数的系数，f()为滤波核函数，所以需满足在N次迭代后，当 k_i∈[0,m]时，低频成分f(k_i)^N≈1，高频成分f(k_i)^N≈0。因此核函数选择式(7)的形式：

f(k)＝(1-λk)(1-μk)  (7)

其中μ是一个新的负尺度因子，绝对值小于1，且μ<-λ。这相当于用正的尺度因子通过式 (1)进行高斯平滑后，进行了另外一个相似的步骤：

x′_i＝x_i+μΔx_i  (8)

由于f(0)＝1，λ+μ＜0，寻找一个临界值k_PB使f(k_PB)＝1，λ、μ满足式(9)：

$k_{\mathrm{PB}}=\frac{1}{\lambda }+\frac{1}{\mu }>0---\left(9\right)$

用式(1)进行高斯平滑时，由于都是加法操作，会导致图形变形，若加入了(8)式的 减法操作，两种计算交替进行，可以克服高斯滤波的该缺陷。这是传统λ/μ方法的原理，它 在高斯滤波基础上修改了核函数，采用加减交替计算的方法，抑制了高斯滤波产生的收缩现 象，但是由于μ<-λ，加运算和减运算的程度不同，所以会影响去噪结果，本发明对λ/μ方法 进行了改进，

当操作数为奇数时进行以下计算：

x′_i＝x_i+λΔx_i，当i为奇数时

x′_i＝x_i+μΔx_i，当i为偶数时

当操作数为偶数时进行如下计算：

x′_i＝x_i+λΔx_i，当i为偶数时

x′_i＝x_i+μΔx_i，当i为奇数时

传统的λ/μ方法每次处理都是方向单一的加权修正，而相邻的点修正方向不一样，这样 多次处理后虽然点与邻域平均值的差异趋小，但相邻两点的差异有可能趋大，造成噪声消除 力度不够。而改进后的方法采用随着处理次数增加而交替使用传统方法，即相邻次数的处理 中，同一个点受到的修正是不同方向的，这样迭代处理，方法本身引入的变形就可以被抵消， 进一步提高去噪效果。

具体实施时，本发明技术方案可采用软件技术实现自动流程运行。本发明实施例所提供 的一种保形LiDAR波形去噪方法，包括如下步骤：

步骤1，对于一维信号x，设包括n个点x₁,…,x_n，标记x_i为一维信号x中第i个点，令操作数 m＝1；

步骤2，进行第m次去噪过程，包括依次对每个点x_i，先通过对邻域计算得出该点的估算 值将x_i原始值与其进行比较和判定，若大于所设定的阈值，则将赋给点x_i， 否则保持原始值；i的取值为2,3,…,n-1，x₁和x_n保持不变，具体实施时，可对x₂,…,x_n-1依次 按此处理，x_i的处理时可以采用x_i-1处理后的结果；

估算值的获取方式为，首先判断当前操作数为奇数还是偶数，

当操作数为奇数时进行以下计算，

x′_i＝x_i+λΔx_i，当i为奇数时

x′_i＝x_i+μΔx_i，当i为偶数时

当操作数为偶数时进行如下计算，

x′_i＝x_i+λΔx_i，当i为偶数时

x′_i＝x_i+μΔx_i，当i为奇数时

其中，变量λ是一个大于0小于1的尺度因子，μ是一个负尺 度因子，且μ<-λ；将计算所得x′_i作为估算值

步骤3，判断当前的操作数m是否达到预设的迭代次数，是则结束流程，得到去噪结果； 否则令m＝m+1，返回步骤2进行迭代，执行下一次去噪过程。

本发明还相应提供一种保形LiDAR波形去噪系统，包括如下模块：

初始化模块，用于对于一维信号x，设包括n个点x₁,…,x_n，标记x_i为一维信号x中第i个点， 令操作数m＝1；

去噪模块，用于进行第m次去噪过程，包括依次对每个点x_i，先通过对邻域计算得出该点 的估算值将x_i原始值与其进行比较和判定，若大于所设定的阈值，则将赋 给点x_i，否则保持原始值；i的取值为2,3,…,n-1，x₁和x_n保持不变；

估算值的获取方式为，首先判断当前操作数为奇数还是偶数，

当操作数为奇数时进行以下计算，

x′_i＝x_i+λΔx_i，当i为奇数时

x′_i＝x_i+μΔx_i，当i为偶数时

当操作数为偶数时进行如下计算，

x′_i＝x_i+λΔx_i，当i为偶数时

x′_i＝x_i+μΔx_i，当i为奇数时

其中，变量λ是一个大于0小于1的尺度因子，μ是一个负尺 度因子，且μ<-λ；将计算所得x′_i作为估算值

迭代判断模块，用于判断当前的操作数m是否达到预设的迭代次数，是则结束工作，得 到去噪结果；否则令m＝m+1，命令去噪模块进行迭代工作，执行下一次去噪过程。

为了验证本发明效果，进行了实验：

1.模拟数据实验

取滤波窗口为3时，分别用均值滤波、高斯滤波、传统λ/μ滤波方法以及改进λ/μ滤波 方法对含有椒盐噪声的模拟LiDAR波形数据进行滤波的结果，其中取k_PB＝0.1，λ＝0.6307，μ ＝-0.6372。滤波窗口取为3时各滤波方法均不能很好地去除噪声，波形中的毛刺依旧很明显， 同时均值滤波和高斯滤波已经开始引起波形的畸变，为达到明显的去噪效果，采用加大滤波 窗口的办法，取滤波窗口为5。加大滤波窗口后，中值滤波、高斯滤波和λ/μ滤波方法对椒 盐噪声的去除有了更明显的效果，但是经均值滤波和高斯滤波去噪后的波形波峰值减小半波 宽增大，发生严重畸变。用传统λ/μ滤波方法可使波形平滑，仅造成了波峰的轻微减小，因 为λ/μ滤波方法在高斯函数的基础上增加了一个抑制波峰变小的因子，可有效克服高斯滤波 导致波形边缘收缩的缺陷，能够在有效去除随机噪声的同时保持波形的形状。从表1中的方 差可以看出改进的λ/μ滤波方法比传统方法去噪效果更好，能够使波形更加接近理想波形。

表1不同方法去噪后波形的方差

滤波方法 均值滤波 高斯滤波 传统λ/μ滤波 改进λ/μ滤波 方差 2.38 6.18 2.07 1.20

2.真实数据实验

此处所采用的真实波形数据是由Leica ALS60系统采集以LAS1.3格式存储的。每个波形 的采样数为128，采样间隔为1纳秒。由上述模拟实验发现，取滤波窗口大小为3时，无论 何种滤波方法都不能有效去除噪声，所以在对真实数据进行实验时，直接取滤波窗口大小为 5。结果中，均值滤波和高斯滤波引起了明显的波形畸变，峰值严重减小。传统和改进的λ/μ 滤波方法的去噪效果从图中看不出明显的区别，进而通过峰值信噪比(PSNR)对两种方法进 行对比，峰值信噪比的定如公式(1)所示，其中max是波形中的最大值，x_i和分别表示第 i个点去噪前和去噪后的值。去噪后波形的峰值信噪比如表2所示，从中可以看出，改进的λ/μ 滤波方法去噪后的波形峰值信噪比更大。

$\mathrm{PSNR}=10\times \lg \frac{{\Sigma }_{i=1}^{128}{\max }^2}{{\Sigma }_{i=1}^{128}{(x_i^{*}-x_i)}^2}---\left(1\right)$

表2各方法对真实数据去噪后的峰值信噪比

滤波方法 均值滤波 高斯滤波 传统λ/μ滤波 改进λ/μ滤波 PSNR 36.24 31.76 40.66 43.95

本文中所描述的具体实施例仅仅是对本发明精神做举例说明。本发明技术领域的技术人 员可以对所描述的具体实施例做各种各样的修改补充或者采用类似的方式替代，但并不会偏 离本发明精神或者超越所附权利要求书所定义的范围。