結合改進DBSCAN和統計濾波的單光子去噪算法

魏 碩，趙楠翔*，李敏樂，胡以華

(1.國防科技大學 電子對抗學院 脈沖功率激光技術國家重點實驗室，合肥 230037；2.國防科技大學 電子對抗學院 先進激光技術安徽省實驗室，合肥 230037)

引 言

近些年,單光子雷達技術逐漸發展成熟[1-3]，作為一種新型激光3維雷達，光子計數激光雷達擁有自己顯著的特點，如回波信號強度低、噪聲大于激光回波、光電探測系統處于光子水平等[4]。不同于傳統雷達探測回波時需要較多的光子形成弱光信號，單光子雷達使用單光子敏感探測器，僅需要檢測到若干光子的返回信息便可得到目標位置信息，使得發射的光子更有效被利用，從而減少激光發射脈沖所需能量，因此單光子探測技術在采集數據方面具有較大優勢，在當下以及未來，光子計數激光雷達在偵察、搜索等方面會發揮越來越重要的作用。

光子計數型激光雷達引入單光子探測器件和時間相關單光子計數 (time-correlated single-photon counting，TCSPC)技術，使系統具有探測極微弱光信號的能力和皮秒量級的時間分辨率[5]，即便在系統設計中采用超窄帶濾波，但由于系統本身的靈敏度極高，探測目標背景仍會產生大量的噪聲，因此，如何有效地去除噪聲是發揮光子計數雷達探測能力的關鍵。目前常用的噪聲去除和數據提取算法主要有殘差比較算法、泊松濾波算法和距離相關檢測算法等[6-8]。

HORAN等人[9]根據掃描數據中信號光子和噪聲光子的分布特征，結合地形設定閾值，從而分離出目標點。HERZFELD等人[10]為模擬的ICESat-2數據分析而開發的計算算法，主要是借助特征向量計算出各項異性的密度中心，再根據地面和樹冠設定密度閾值，達到識別地物的目的。LI等人[11]對初步去噪后的點云使用改進后的統計濾波方法，設定閾值進行去噪，但參量的設定比較依賴經驗，對結果也有著較大的影響。LI等人[12]依靠計算點云密度進行去噪，首先進行圖像分割進行粗去噪，再定義一種線性徑向基函數進行精去噪，得到結果與設定的閾值比較，得出較好的結果。XU等人[13]在對單光子激光數據特點精確分析的基礎上，提出一種基于地形相關和最小二乘曲線擬合的數據處理方法，獲得較為精確的結果。上述算法要求的參量較多，且參量數值選取有較高的要求，作者在研究單光子數據空間分布特征的基礎上，提出一種基于密度、參量自適應較強的去噪算法，并通過美國國家航空航天局(National Aeronautics and Space Administration ，NASA)提供的多波束試驗激光雷達(multiple altimeter beam experimental lidar,MABEL)實際飛行數據[14-15]進行驗證，得到較好的結果。

1 光子計數激光雷達系統原理

光子計數激光雷達[16]的系統框圖如圖1所示。

Fig.1 Photon counting lidar system

如圖1所示，光子計數激光雷達的工作原理是：脈沖激光器發射的脈沖經分光片后分成2束，一束充當計時脈沖照射到光電探測器，觸發計時信號；其余部分照射到物體，經反射后產生回波被單光子探測器接收，并在計時電路中產生停止信號。TCSPC模塊計算2次信號的時間差，得到本次脈沖的來回飛行時間。隨著探測時間和探測次數的增加，得到掃描光子隨時間分布的數據，傳遞給后續模塊通過時間獲取算法得到具體目標距離數據及成像。

由于系統本身的靈敏度很高，在探測過程中目標背景會產生大量的噪聲，激光雷達的噪聲大致分為3類：探測器噪聲、放大器噪聲、背景的輻射噪聲[17]。隨著制造工藝和致冷方式的進步，探測器噪聲和放大器噪聲得到有效抑制，為抑制背景噪聲，除了在接收光學系統添加窄帶濾波片外，還會采用距離門技術。

距離門技術是指給探測器添加一個時間區間，只有在這個時間區間內返回的回波光子才可以被探測到，在這個時間段外的回波光子不會被探測到。雖然如此，但大量的背景噪聲仍會被引入到探測數據中，這就需要更近一步的去噪[18-19]。

2 本文中的算法

2.1 算法相關概念

2.1.1 點密度xi常規的點云密度定義Di是指在規定的半徑eps內包含點云的數目，如下式所示：

Di=number(|pi-pj|

(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個點云的最小距離。

xi=max{x1,x2,…,xk}

(2)

2.1.2 密度閾值xm分割閾值xm根據點密度xi求出，其大小為整體點密度的均值，計算下式所：

(3)

由光子計數激光雷達的特性以及對實際數據集的研究可知，噪聲點云的數量遠遠大于信號點云，因此分割閾值的范圍在0.8xm～xm之間均可以滿足實驗要求，經多次實驗可知，分割閾值的取值為0.9xm時結果最佳。

2.1.3 改進的DBSCAN算法 DBSCAN算法根據人為輸入參量eps以及minpts(DBSCAN算法中的固定變量符號，表示eps距離內包含的個數)衡量點云密度，遍歷所有點云，將點云劃分成一個個密度不等的簇。由于輸入參量有較大的主觀性，如果和數據集密度特征不匹配容易造成過分割或欠分割。

為了彌補原DBSCAN算法的缺點，本文中將點密度xi用于衡量單個點云密度，由于點密度xi和密度閾值xm均由數據集本身密度特征自動計算，因此在核心點“生長”為簇的過程中，避免了人工帶來的誤差。

改進后的聚類算法由于具有密度參量自適應的特性，可以對絕大多數的點云集進行合理分割。

2.2 算法詳細步驟

本文中去噪方法可以分為2步：粗去噪和精去噪。首先對數據進行粗去噪，去除大部分明顯的噪聲點，可以極大地簡化后續算法的計算量；然后根據剩余點云的密度特征和空間分布特征，執行精去噪，獲得目標點云數據。算法流程如圖2所示。

Fig.2 The flow chart of noise filter method

2.2.1 粗去噪 (1)坐標轉換。實際實驗數據是大地坐標系(L,M,H)，即由經度L、維度M和高程H組成，利用坐標轉換[20-21]，將大地坐標系數據轉化成空間直角坐標系(x,y,z)數據；(2)建立空間拓撲關系。對直角坐標系下的點云數據通過k維樹建立空間拓撲關系，由近鄰索引計算點云之間密度；(3)設定閾值去除明顯噪聲。選取合適的閾值進行點云分割(本文中選取0.9倍密度均值作為閾值)，濾除明顯的離散噪聲點。

di=

(4)

(5)

(6)

3 參量自適應

本文中的算法，需要設定的參量只有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

此外還借助2個參考標準衡量k值對去噪結果的影響。一個是濾波精度N，其公式如下：

N=Nc/Nt

(7)

式中，Nc表示正確區分的信號點數目，Nt表示信號點總數。濾波精度N主要描述了點云正確識別率，但沒有考慮噪聲的影響。基于此，本文中還采用交并比U來描述信號點與噪聲點的關系，公式如下：

U=(A∩B)/(A∪B)

(8)

式中,A表示實際信號點云數據集，B表示去噪后點云數據集，(A∩B)表示的是去噪后保留的點云數，(A∪B)表示的是實際信號點云數據和去噪后保留的噪聲數據之和。一般來說，在計算機檢測任務中，U>0.5的結果就是可以接受的，U越大，則點云圖的質量也就越好。具體數據如表2所示。

Table 2 The relationship between k value and N,U

從表1可知，k值和去噪后的信號點和噪聲點數目均呈正相關，這和上面對k值的分析相符。由表2可知，隨著k值的增長，N的增長率逐漸接近于0，而U卻逐漸降低，表明噪聲點的增長率大于信號點的增長率，繼續增大k值得不償失。總體上來看，U保持在75%～80%左右，表明k值的適當選取對點云圖質量影響不是很大，符合本文中算法參量自適應的要求。

4 實驗結果與討論

4.1 實驗結果

由第3節中分析可知，本文中的算法參量具有參量自適應強的特點，k值對最終結果的影響并不是很大，因此令k=23，首先將數據由大地坐標系轉換到空間直角坐標系，并將坐標系原點平移到數據最小處，然后應用本文中的算法進行去噪。去噪后，剩余噪聲點為6290個，保留的目標點為41733個，濾波精度N=90.27%，交并比U=79.5%。

圖3是坐標轉換且進行平移后的原始數據。其中黑色是噪聲，細線是山脈。從圖中可以看出,噪聲點的數量遠大于信號點數量。圖4是粗去噪的結果。從圖中可以看出，大量的噪聲點已經被去掉，可以近似看出山脈的輪廓。圖5和圖6是精去噪的結果圖。其中圖5是利用改進DBSCAN算法進行二次去噪后的結果，圖中極大部分明顯的團狀噪聲點已經被去除。圖6是應用統計濾波處理后的最終去噪結果，目標主體附近的離散噪聲點云也有了明顯的去除。圖7是部分放大后的場景圖。從中可以看出，目標細節有了更多的保留，山坡上的樹木也可以分辨。

Fig.3 3-D scene graph

Fig.4 Rough denoising result graph

Fig.5 Second denoising result graph

Fig.6 The final result graph

Fig.7 Partially enlarged view

4.2 算法對比

半徑濾波算法也是常用的濾波算法，其原理是假定數據集中每一個點在給定的半徑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

結合圖8和表3可以發現，在不同的參量下，半徑濾波的結果起伏非常大。半徑濾波的準確性建立在較為苛刻的先驗知識上，需要人為地設定參量，不具備普遍適用性，最好的參量情況下，其結果達到78.39%。

參考文獻[12]中提出的去噪算法也采用濾波精度作為衡量算法的指標，該文獻中濾波精度分別達到90.26%以及96.11%。本文中算法在輸入參量k>20后，其濾波精度均達到90%左右，點云去噪結果較好，且可選參量范圍較廣，突出本文中算法參量自適應能力強的特點。

Fig.8 Filtering results under different r valuesa—filtering results under r=9 b—filtering results under r=12

5 結 論

DBSCAN算法分割結果的好壞取決于參量設定，需要eps,minpts和數據集密度特征相匹配，因此，DBSCAN算法在處理不同數據集時需要多次實驗來確定最佳參量，而本文中改進后的DBSCAN算法僅需要設定臨近點數目，相應的距離參量可以由k維樹根據數據集密度特征自動得出，因此改進的DBSCAN算法在適應性上大大提高。其次，本文中對改進后DBSCAN算法分割后的結果進行處理，去除了目標附近大部分團狀噪聲點云，再根據剩余點云特征，設定統計濾波法的閾值對點云進行最終處理，得到結果。

半徑濾波算法是應用極為廣泛的密度濾波算法，本文中利用半徑濾波算法處理數據集，將得到結果與本文中算法相對比，可以看出,本文中算法幾乎不需要先驗知識，在參量自適應性上遠遠優于半徑濾波算法。此外,本文中的去噪算法在參量適中的情況下,U普遍在75%～80%(半徑濾波算法的U在20%～79%),濾波精度在85%～93%，去噪結果也好于半徑濾波算法。兩種算法對比，更加突出本文中算法的優越性。

但本文中的算法也存在一定的缺點，對于距離目標較近且密度極大的團狀點云無法有效的去除，只能盡量減少其大小與數量，此外，如果數據集很大的情況下，k維樹的搜索速度也逐漸下降。由圖5粗去噪結果圖可以清晰地看出掃描線結構，在粗去噪的基礎上，結合掃描線間距以及同一掃描線上點間距進行精確去噪，這是作者下一步進行改進的方向。