-
常规的点云密度定义Di是指在规定的半径eps内包含点云的数目,如下式所示:
$ D_{i}=\operatorname{number}\left(\left|p_{i}-p_{j}\right|<\text { eps }\right) $
(1) 式中, pi是中心点坐标,pj是临近点坐标,i, j=1, 2, 3…; eps是指基于密度的空间聚类噪声应用(density-based spatial clustering of applications with noise,DBSCAN)算法中的半径参量。
与常规的点云密度定义相反,本文中规定点密度xi表示单个点云的密度,如(2)式所示。由k维树可以得到k个临近点到中心点的距离数组{x1, x2, …, xk},由于k值是输入参量,因此点密度xi表示包含k个点云的最小距离。
$ x_{i}=\max \left\{x_{1}, x_{2}, \cdots, x_{k}\right\} $
(2) -
分割阈值xm根据点密度xi求出,其大小为整体点密度的均值,计算下式所:
$ x_{\mathrm{m}}=\frac{x_{1}+x_{2}+\cdots+x_{n}}{n} $
(3) 由光子计数激光雷达的特性以及对实际数据集的研究可知,噪声点云的数量远远大于信号点云,因此分割阈值的范围在0.8xm~xm之间均可以满足实验要求,经多次实验可知,分割阈值的取值为0.9xm时结果最佳。
-
DBSCAN算法根据人为输入参量eps以及minpts(DBSCAN算法中的固定变量符号,表示eps距离内包含的个数) 衡量点云密度,遍历所有点云,将点云划分成一个个密度不等的簇。由于输入参量有较大的主观性,如果和数据集密度特征不匹配容易造成过分割或欠分割。
为了弥补原DBSCAN算法的缺点,本文中将点密度xi用于衡量单个点云密度,由于点密度xi和密度阈值xm均由数据集本身密度特征自动计算,因此在核心点“生长”为簇的过程中,避免了人工带来的误差。
改进后的聚类算法由于具有密度参量自适应的特性,可以对绝大多数的点云集进行合理分割。
-
本文中去噪方法可以分为2步:粗去噪和精去噪。首先对数据进行粗去噪,去除大部分明显的噪声点,可以极大地简化后续算法的计算量;然后根据剩余点云的密度特征和空间分布特征,执行精去噪,获得目标点云数据。算法流程如图 2所示。
-
(1) 坐标转换。实际实验数据是大地坐标系(L, M, H),即由经度L、维度M和高程H组成,利用坐标转换[20-21],将大地坐标系数据转化成空间直角坐标系(x, y, z)数据;(2)建立空间拓扑关系。对直角坐标系下的点云数据通过k维树建立空间拓扑关系,由近邻索引计算点云之间密度;(3)设定阈值去除明显噪声。选取合适的阈值进行点云分割(本文中选取0.9倍密度均值作为阈值),滤除明显的离散噪声点。
-
(1) 计算剩余数据的点密度xi。对剩下的数据建立k维树,由最近邻索引得到距离数组{x1, x2, …, xk},求出点密度xi;(2)根据点密度进行点云分簇。随机选取一个点,当该点点密度满足阈值xm时进行区域增长,逐步遍历其周围的点,将满足条件的点归为一簇,重复此过程,直至所有点云均被处理;(3)去除团状噪声簇。数据分簇后,检查每个点云簇包含的点云数量,对于点云数小于一定值的点云簇视为噪声簇去除,具体分析如下:由于本文中选取的研究对象是山脉,在噪声影响的下,目标也可能被分成几段,所以经过试验选取阈值为0.5T,既可以达到去噪目的也尽可能保留信号点,其中,阈值T=数据集点云数÷聚类后簇数,查看每个簇点云数目,如果小于0.5T,则作为密集噪声簇去除; (4)去除目标很近的离散噪声。遍历点云得到临近的k个点, 并求出距离均值构成1维数组D,然后计算数组D的均值di和标准差σ,设定阈值滤除离散点。
定义di为点pi(xi, yi, zi)到其k个临近点pij(xij, yij, zij)的平均距离,di为di的均值,σ为di的标准差,j=1, 2, …, k。
$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{d_i} = \\ \frac{{\sum\limits_{j = 1}^k {\sqrt {{{\left( {{x_{ij}} - {x_i}} \right)}^2} + {{\left( {{y_{ij}} - {y_i}} \right)}^2} + {{\left( {{z_{ij}} - {z_i}} \right)}^2}} } }}{k} \end{array} $
(4) $ {{\bar d}_i} = \frac{{\sum\limits_{i = 1}^n {{d_i}} }}{n} $
(5) $ \sigma = \sqrt {\frac{1}{n}\sum\limits_{i = 1}^n {{{\left( {{d_i} - {{\bar d}_i}} \right)}^2}} } $
(6) -
本文中的算法,需要设定的参量只有k,k表示选取临近点的数目,k值大小的意义在于点密度的衡量标准。k值很小, 话点密度由最近的较少点决定,增加了特殊值出现的情况,如某个信号点恰好周围最近的几个点较稀疏,或者某个噪声点最近的几个点距离很近等;k值很大, 会使信号点与噪声点的密度分割线较为模糊,虽然可以额外保留分布特殊的信号点,但也容易保留较多的噪声点。
本文中的数据来源是美国NASA提供的MABEL实际飞行数据,目标区域是位于美国加州Sierras-Forest地区,采用2010年12月飞行试验的532nm波段第50通道的观测数据进行实验。坐标转化后的数据总共有169511个点,其中噪声数据点有123280个,目标点云有46231个。
本文中选取不同k值代入算法进行运算,去噪后的点云数据如表 1所示。
Table 1. Point cloud data after denoising
k 8 12 16 25 38 50 70 100 noise point 7309 5592 5309 7205 8064 9411 9859 11044 signal point 40008 40575 40867 42741 42231 42674 42283 42822 此外还借助2个参考标准衡量k值对去噪结果的影响。一个是滤波精度N,其公式如下:
$ N=N_{\mathrm{c}} / N_{\mathrm{t}} $
(7) 式中,Nc表示正确区分的信号点数目,Nt表示信号点总数。滤波精度N主要描述了点云正确识别率,但没有考虑噪声的影响。基于此,本文中还采用交并比U来描述信号点与噪声点的关系,公式如下:
$ U=(A \cap B) /(A \cup B) $
(8) 式中, A表示实际信号点云数据集,B表示去噪后点云数据集,(A∩B)表示的是去噪后保留的点云数,(A∪B)表示的是实际信号点云数据和去噪后保留的噪声数据之和。一般来说,在计算机检测任务中,U>0.5的结果就是可以接受的,U越大,则点云图的质量也就越好。具体数据如表 2所示。
Table 2. The relationship between k value and N, U
k 8 12 16 25 38 50 70 100 N/% 86.6 87.8 88.4 92.5 91.2 92.3 91.5 92.6 U/% 74.7 78.3 79.3 80.0 77.8 76.7 75.4 74.8 从表 1可知,k值和去噪后的信号点和噪声点数目均呈正相关,这和上面对k值的分析相符。由表 2可知,随着k值的增长,N的增长率逐渐接近于0,而U却逐渐降低,表明噪声点的增长率大于信号点的增长率,继续增大k值得不偿失。总体上来看,U保持在75%~80%左右,表明k值的适当选取对点云图质量影响不是很大,符合本文中算法参量自适应的要求。
-
由第3节中分析可知,本文中的算法参量具有参量自适应强的特点,k值对最终结果的影响并不是很大,因此令k=23,首先将数据由大地坐标系转换到空间直角坐标系,并将坐标系原点平移到数据最小处,然后应用本文中的算法进行去噪。去噪后,剩余噪声点为6290个,保留的目标点为41733个,滤波精度N=90.27%,交并比U=79.5%。
图 3是坐标转换且进行平移后的原始数据。其中黑色是噪声,细线是山脉。从图中可以看出, 噪声点的数量远大于信号点数量。图 4是粗去噪的结果。从图中可以看出,大量的噪声点已经被去掉,可以近似看出山脉的轮廓。图 5和图 6是精去噪的结果图。其中图 5是利用改进DBSCAN算法进行二次去噪后的结果,图中极大部分明显的团状噪声点已经被去除。图 6是应用统计滤波处理后的最终去噪结果,目标主体附近的离散噪声点云也有了明显的去除。图 7是部分放大后的场景图。从中可以看出,目标细节有了更多的保留,山坡上的树木也可以分辨。
-
半径滤波算法也是常用的滤波算法,其原理是假定数据集中每一个点在给定的半径r内至少存在k个点,符合假定条件的作为信号点保存下来,不符合条件的作为噪声点除去,其中半径r以及数值k由人工指定。为了优化结果,事先选取数值k=30,半径r在密度阈值附近多个范围内取值,观察不同参量下的半径滤波结果,并与本文中去噪算法结果进行比较。
在k=30时,计算出的r=21。由于r值越大,保留的噪声点就越多,因此对于参量r主要选取小于21的值。令r为3, 6, 9, 12, 15, 18, 21和24,去噪后数据如表 3所示,去噪效果较好的结果如图 8所示。
Table 3. Radius filtering data
r 3 6 9 12 15 18 21 24 signal point 9542 32291 40673 44049 45481 45990 46165 46198 noise point 29 1522 5652 10425 19011 31771 44969 58410 U/% 20.62 67.62 78.39 77.74 69.71 58.96 50.61 44.14 结合图 8和表 3可以发现,在不同的参量下,半径滤波的结果起伏非常大。半径滤波的准确性建立在较为苛刻的先验知识上,需要人为地设定参量,不具备普遍适用性,最好的参量情况下,其结果达到78.39%。
参考文献[12]中提出的去噪算法也采用滤波精度作为衡量算法的指标,该文献中滤波精度分别达到90.26%以及96.11%。本文中算法在输入参量k>20后,其滤波精度均达到90%左右,点云去噪结果较好,且可选参量范围较广,突出本文中算法参量自适应能力强的特点。
结合改进DBSCAN和统计滤波的单光子去噪算法
Single photon denoising algorithm combined with improved DBSCAN and statistical filtering
-
摘要: 为了解决光子计数激光雷达探测数据中噪声点云过多的问题, 采用结合基于密度的噪声空间聚类应用算法(DBSCAN)和统计滤波算法的单光子点云去噪方法, 以美国国家航空航天局提供的多波束试验激光雷达实际飞行数据为实验数据, 通过k维树求取点云密度进行粗去噪, 然后运用改进DBSCAN算法和统计滤波算法进行精去噪, 进行了理论分析和实验验证。结果表明, 实验区目标点云识别率在85%以上, 性能优于经典的半径滤波算法。这一结果对于光子数据去噪是有帮助的。Abstract: In order to solve the problem of excessive noise point clouds in the photon counting lidar detection data, a single photon point cloud denoising method based on a combination of improved density-based spatial clustering of applications with noise (DBSCAN) algorithm and statistical filtering algorithm was adopted. The actual flight data of multiple altimeter beam experimental lidar provided by National Aeronautics and Space Administration was experimental data. First, the point cloud density was obtained through the k-dimensional tree for rough denoising, and then the improved DBSCAN algorithm and statistical filtering algorithm were used for fine denoising. The theoretical analysis and experimental verification has achieved good results. The results show that the target point cloud recognition rate in the experimental area is above 85%, and the performance is better than the classic radius filtering algorithm. This result is helpful for photon data denoising.
-
Table 1. Point cloud data after denoising
k 8 12 16 25 38 50 70 100 noise point 7309 5592 5309 7205 8064 9411 9859 11044 signal point 40008 40575 40867 42741 42231 42674 42283 42822 Table 2. The relationship between k value and N, U
k 8 12 16 25 38 50 70 100 N/% 86.6 87.8 88.4 92.5 91.2 92.3 91.5 92.6 U/% 74.7 78.3 79.3 80.0 77.8 76.7 75.4 74.8 Table 3. Radius filtering data
r 3 6 9 12 15 18 21 24 signal point 9542 32291 40673 44049 45481 45990 46165 46198 noise point 29 1522 5652 10425 19011 31771 44969 58410 U/% 20.62 67.62 78.39 77.74 69.71 58.96 50.61 44.14 -
[1] XU Q, XIE X M, ZHANG W, et al. The progress of semiconductor quantum dot based quantum emitter [J]. Laser Technology, 2020, 44(5): 575-586(in Chinese). [2] TANG X M, LI G Y. Development and prospect of laser altimetry satellite [J]. International Space, 2017(11): 13-18(in Chinese). [3] LUO H J, ZHOU R L, ZHANG Y T. Theoretical analysis of detection performance and range accuracy of photon ladar[J]. Laser Technology, 2014, 38(3): 411-416(in Chinese). [4] SHU R, HUANG G H, HOU L B, et al. Multi-channel photon counting three-dimensional imaging laser radar system using fiber array coupled Geiger-mode avalanche photodiode[J]. Proceedings of the SPIE, 2012, 8542: 1-10. [5] FANG J, SHE Ch, LIU J P. A denoising method based on photon counting lidar[J]. Ship Electronic Warfare, 2019, 42(4): 10-15(in Chinese). [6] OH M S, KONG H J, KIM T H, et al. Multihit mode direct-detection laser radar system using a Geiger-mode avalanche photodiode[J]. Review of Scientific Instruments, 2010, 81(3): 1-7. [7] FOUCHE D G. Detection and false-alarm probabilities for laser radars that use Geiger-mode detectors[J]. Applied Optics, 2003, 42(27): 5388-5398. doi: 10.1364/AO.42.005388 [8] MILSTEIN A B, JIANG L A, LUU J X. Acquisition algorithm for direct-detection ladars with Geiger-mode avalanche photodiodes[J]. Applied Optics, 2008, 47(2): 296-311. doi: 10.1364/AO.47.000296 [9] HORAN K H, KEREKES J P. An automated statistical analysis approach to noise reduction for photon-counting lidar systems [C]//IEEE International Geoscience and Remote Sensing Symposium. NewYork, USA: IEEE, 2013: 4336-4339. [10] HERZFELD U C, MCDONALD B W, WALLIN B F, et al. Algorithm for detection of ground and canopy cover in micropulse photon-counting lidar altimeter data in preparation for the ICESat-2 mission [J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(4) : 2109-2125. doi: 10.1109/TGRS.2013.2258350 [11] LI K, ZHANG Y Sh, TONG X Ch, et al. Research on de-noising and filtering algorithm of single photon lidar data[J]. Navigation and Control, 2020, 19(1): 67-76(in Chinese). [12] LI M, GUO Y, YANG G, et al. A push-broom photon counting lidar point cloud filtering algorithm and its verification[J]. Science Technology and Engineering, 2017, 17(9): 53-58 (in Chinese). [13] XU Y T, LI G Y, QIU Ch X, et al. Single-photon laser data processing technology based on terrain correlation and least square curve fitting[J]. Infrared and Laser Engineering, 2019, 48(12): 148-157(in Chinese). [14] MCGILL M, MARKUS T, SCOTT V S, et al. The multiple altimeter beam experimental lidar (MABEL): An airborne simulator for the ICESat-2 mission[J]. Journal of Atmospheric and Oceanic Technology, 2013, 30(2): 345-352. doi: 10.1175/JTECH-D-12-00076.1 [15] BRUNT K M, NEUMANN T A, AMUNDSON J M, et al. MABEL photon-counting laser altimetry data in Alaska for ICESat-2 simulations and development[J]. Cryosphere, 2016, 10(4): 1707-1719. doi: 10.5194/tc-10-1707-2016 [16] YU F L. Research on ladar 3-D point cloud data processing methods based on single photon detection [D]. Harbin: Harbin University of Technology, 2016: 57-63(in Chinese). [17] ZHOU H, LI S, WANG L X, et al. Influence of noise on range error for satellite laser altimeter [J]. Infrared and Laser Engineering, 2015, 449(8): 2256 - 2261(in Chinese). [18] LIU Y, SUN Sh Y. Laser point cloud denoising based on principal component analysis and surface fitting[J]. Laser Technology, 2020, 44(4): 497-502(in Chinese). [19] WANG T, SHEN Y H, YAO J Q. Research on laser radar echo signal denoising based on wavelet threshold method[J]. Laser Technology, 2019, 43(1): 63-68(in Chinese). [20] LIU Zh P, YU Q Y, CHA J F. Rapid transformation from spatial rectangular coordinates to two common coordinates[J]. Science of Surveying and Mapping, 2015, 40(3): 8-11(in Chinese). [21] SANG J. Non-iteration method for inversion of space geodetic rectangular coordinates and geodetic coordinates[J]. Bulletin of Surveying and Mapping, 2000(5): 37-39(in Chinese).