劉翔，張立華，戴澤源，陳秋，周寅飛

（1 海軍大連艦艇學院 軍事海洋與測繪系，遼寧 大連 116018）（2 海軍大連艦艇學院海洋測繪工程軍隊重點實驗室，遼寧 大連 116018）（3 91001部隊，北京 100071）

0 引言

2018年9月15日發(fā)射的ICESat-2（The Ice，Cloud and Land Elevation Satellite-2）是美國第二代冰、云和陸地高程衛(wèi)星，用以執(zhí)行冰蓋高程測量和變化監(jiān)測、陸地高程測量、全球植被高度測量以及監(jiān)測云和氣溶膠等任務(wù)［1-3］。相較于ICESat所搭載的全波形地學激光測高系統(tǒng)（Geoscience Laser Altimeter System，GLAS），ICESat-2搭載的先進地形激光高度計系統(tǒng)（Advanced Topographic Laser Altimeter System，ATLAS）為微脈沖光子計數(shù)測高雷達系統(tǒng)，具有高靈敏性、高重頻的特性，可以獲得光斑更小、密度更高的光子點云數(shù)據(jù)，增強了對全球三維信息的獲取能力［2］。同時這也是全球范圍內(nèi)第一次將該項技術(shù)運用于星載平臺之上［1］，代表了未來星載激光測高重要發(fā)展趨勢。但由于ICESat-2光子計數(shù)激光雷達發(fā)射和探測的激光脈沖都為弱信號，難以區(qū)分目標物體表面返回脈沖、大氣散射、太陽輻射及儀器自身噪聲［4-5］，導致其回波數(shù)據(jù)中摻雜了大量的背景噪聲。在強噪聲背景下，較低的信噪比使得信號光子與噪聲光子難以區(qū)分，如何從背景噪聲中提取信號光子已成為ICESat-2數(shù)據(jù)應(yīng)用的一項關(guān)鍵而基礎(chǔ)的工作。

針對單光子激光測高雷達數(shù)據(jù)的去噪，國內(nèi)外學者進行了諸多研究。文獻［5-6］通過統(tǒng)計光子與其鄰近光子局部距離進行去噪處理，但鄰近光子數(shù)量為經(jīng)驗設(shè)定值；文獻［7-11］基于光子濾波核內(nèi)光子密度值的差異進行去噪，但濾波核范圍參數(shù)需人為進行設(shè)定；文獻［12］以DBSCAN聚類方法為基礎(chǔ)，將所用圓形濾波核形狀改為水平橢圓進行信號光子聚類，但該聚類方法存在對輸入?yún)?shù)敏感的問題；文獻［13］則以O(shè)PTICS方法進行聚類，降低了對輸入?yún)?shù)的敏感度，但所采用橢圓濾波核半徑仍需通過試錯法來得到。以上去噪方法存在參數(shù)人為設(shè)定較為繁雜或?qū)?shù)敏感的問題。盡管文獻［14-15］將光子分布光柵化為圖像，以邊緣提取的圖像處理方式進行去噪，避免了參數(shù)設(shè)定過程，但光柵化必然會造成光子點云信息不可逆的損失，影響數(shù)據(jù)的后續(xù)應(yīng)用。

為了解決單光子激光測高雷達數(shù)據(jù)去噪對輸入?yún)?shù)的依賴性問題，ZHANG Guoping等［16］提出了一種基于四叉樹的無參數(shù)輸入去噪方法，該方法無需人為經(jīng)驗設(shè)定的光子在高程和沿軌距離方向上的范圍以及周圍光子數(shù)量，只需利用四叉樹將每個光子合理分隔即可通過分隔層值表征光子密度，從而可做到在去噪時無需要參數(shù)輸入，在ICESat-2的模擬數(shù)據(jù)MATLAS去噪實驗中取得了良好的效果。但對于強噪聲背景下的實測ICESat-2數(shù)據(jù)，四叉樹去噪方法處理局部稀疏、相距較近的噪聲光子時表征密度不合理，會將較多噪聲光子誤識別為信號光子，從而影響去噪效果。

為此，本文針對四叉樹去噪方法的缺陷，首先以剪枝四叉樹的層值表征光子密度，避免強噪聲背景下局部稀疏、相距較近的噪聲光子被四叉樹分隔層次較多而被誤表征為密度過大，自適應(yīng)求取信號光子的密度閾值，完成一級去噪；然后通過箱線圖進一步去除少量未能被剪枝四叉樹識別局部較密離群噪聲光子，完成二級去噪。

1 數(shù)據(jù)與研究方法

1.1 研究區(qū)域概況

本文研究區(qū)域及實驗所用ICESat-2數(shù)據(jù)分布情況如圖1所示，研究區(qū)域1位于美國的北達科他州，該區(qū)域為平原地帶，地勢較為平坦，植被覆蓋稀疏；研究區(qū)域2位于美國的加利福尼亞州，該區(qū)域靠近內(nèi)華達山脈，地勢起伏較大，中等植被覆蓋。兩個研究區(qū)域地勢起伏差異較大、地表覆蓋類型不同，具有較強的代表性。

圖1 研究區(qū)域位置及數(shù)據(jù)分布示意圖Fig.1 Location of study area and distribution of data

1.2 數(shù)據(jù)來源

1.2.1 ICESat-2/ATLAS數(shù)據(jù)

2019年5 月美國國家冰雪數(shù)據(jù)中心NSIDC（National Snow and Ice Data Center）開始向全球用戶公開發(fā)布ICESat-2/ATLAS數(shù)據(jù)［4］。本文使用其中的Level-2級數(shù)據(jù)產(chǎn)品ATL03作為本次去噪的實驗數(shù)據(jù)，ATL03為全球定位光子數(shù)據(jù)，它記錄了光子精確的經(jīng)緯度和高程信息，是生成諸如海冰高程（ATL07）［17］、陸地植被高程（ATL08）［18］等其他產(chǎn)品的數(shù)據(jù)基礎(chǔ)［19］。

1.2.2 機載激光雷達數(shù)據(jù)

為驗證本文方法去噪效果，實驗選擇美國國家生態(tài)觀測站網(wǎng)絡(luò)NEON（National Ecological Observatory Network）發(fā)布的機載激光雷達高程數(shù)據(jù)產(chǎn)品作為驗證數(shù)據(jù)，包括數(shù)字表面模型（Digital Surface Model，DSM）和數(shù)字地面模型（Digital Terrain Model，DTM），二者分辨率皆為1 m。同時為保證驗證的準確性，機載激光雷達數(shù)據(jù)觀測時間與ICESat-2數(shù)據(jù)觀測時間均為同年同月，以減少因時相變化所造成的不同數(shù)據(jù)產(chǎn)品之間的高程誤差。

1.3 研究方法

1.3.1 一級去噪

在ICESat-2光子點云中，信號光子與噪聲光子空間分布差異往往較大，具體表現(xiàn)為信號光子較噪聲光子分布更加密集［20-21］，這種密度的差異是當前ICESat-2點云去噪的基本依據(jù)。因此，選取合適的密度閾值是正確分類信號光子和噪聲光子的關(guān)鍵。本文根據(jù)信號光子和噪聲光子密度不同的特點，挖掘Otsu法在自適應(yīng)求取灰度特性有明顯差別的兩類目標分割閾值上的優(yōu)勢，通過剪枝四叉樹的層值來定量表示光子密度，以密度替代灰度，求取信號光子和噪聲光子密度最大類間方差作為分類閾值，對兩類光子進行分類。

1.3.1.1 四叉樹層值表征光子密度

四叉樹層值表征光子密度是利用四叉樹方法對光子分隔，即在高程和沿軌距離二維空間D上按照光子的位置將空間進行四象限遞歸劃分，依照光子最終所處四叉樹象限葉節(jié)點的層值來定量光子的的密度。分隔過程為首先根據(jù)ICESat-2光子數(shù)據(jù)的范圍給定初始根節(jié)點空間D0，記錄D0中光子個數(shù)為sum(D0)，當sum(D0)超過設(shè)定的空間光子最大容納值maxlimit，則將當前空間劃分為右上（Ⅰ）、左上（Ⅱ）、左下（Ⅲ）、右下（Ⅳ）四個相同的子象限空間，計算各子空間內(nèi)光子數(shù)，對超過maxlimit的子空間再次進行劃分，遞歸執(zhí)行以上劃分，直至各子空間內(nèi)光子數(shù)在最大容納值之內(nèi)，此時光子處于各葉節(jié)點所對應(yīng)的劃分空間中。空間劃分準則可歸納為

式中，D′k-1為第k-1層中可進行劃分的象限集合；Dk-1為第k-1層象限集合；X為Dk-1中某個象限空間；sum(X)為象限空間X內(nèi)光子數(shù)；為保證去噪精度，通常需對每個光子進行分隔，因此maxlimit一般設(shè)為1。

圖2為某區(qū)域利用四叉樹分隔陰影部分光子的示意圖，圖3為相應(yīng)的四叉樹樹狀結(jié)構(gòu)圖，從中可以看出由于噪聲光子分布較信號光子稀疏，因此在分隔時只需要較少的劃分次數(shù)即可被分隔開，處于的葉節(jié)點象限空間層值k較小，而信號光子因分布密集所處葉節(jié)點象限空間層值k較高，層值k與密度存在著明顯的正相關(guān)關(guān)系，因此采用層值來對信號和噪聲光子密度進行表征，將層值作為光子的密度標簽。

圖3 四叉樹分隔光子樹狀結(jié)構(gòu)圖Fig.3 Tree structure of photons separated by quadtree

1.3.1.2 剪枝四叉樹層值表征光子密度

由于噪聲光子的分布呈現(xiàn)隨機性，在利用四叉樹分隔ICESat-2光子時，會出現(xiàn)局部稀疏情形下某些噪聲光子相距較近而需要較多的劃分次數(shù)方可分隔的現(xiàn)象，如圖4紅框標識的象限所示。此時利用四叉樹方法表征光子密度時，這類噪聲光子由于所在葉節(jié)點象限層值較高，接近信號光子的層值，如圖5所示，導致在去噪過程中被誤分為信號光子。特別是在強噪聲背景下，信噪比較低，光子點云中會大量散布此類噪聲光子，進而導致四叉樹方法難以達到較好的去噪效果。

圖5 四叉樹誤分隔光子樹狀結(jié)構(gòu)圖Fig.5 Tree structure of photons incorrectly separated by quadtree

進一步分析光子密度的表征過程，光子層值需要反映的是光子與周圍光子而不是與最鄰近光子之間的疏密程度，在圖4紅框標識的象限中，內(nèi)部光子雖然距離較近，但是象限進行一次四等劃分后并不能將其分隔，從其與周圍光子分布的稀疏程度來看其密度仍然較小，所以此時象限空間所在層對應(yīng)的層值就代表了其內(nèi)部光子的密度。基于此，為減少噪聲光子的誤識別，本文引入數(shù)據(jù)搜索過程中的“剪枝”思想改進原有四叉樹生成過程，并采用剪枝四叉樹層值來表征光子密度。數(shù)據(jù)搜索算法中的“剪枝”是通過一定的策略對樹結(jié)構(gòu)進行優(yōu)化，刪除部分無效節(jié)點，提升算法效率。而本文所述的“剪枝”過程，則是先對將要進行劃分的某象限進行分析判斷，若劃分一次后沒有光子被分隔開，則此劃分不進行，同時此象限的劃分相應(yīng)終止，以此來剪掉由于不合理的劃分而生成的枝狀結(jié)構(gòu)。劃分準則為

圖4 四叉樹誤分隔光子示意圖Fig.4 Photons incorrectly separated by quadtree

式中，nsum[split(X)]≥1表示為對象限空間X四等劃分后，子象限空間內(nèi)光子數(shù)不小于1的空間個數(shù)，即存在光子的子象限空間個數(shù)。

對四叉樹劃分改進后，劃分準則不再依據(jù)maxlimit，而是顧及了象限空間光子的分布情況。按照劃分準則，對象限空間X的劃分存在以下兩種情況：

1）sum(X)≤1，即象限空間X不存在或只存在一個光子，對此象限不進行劃分，此象限劃分終止；

2）sum(X)≥2，即象限空間X存在多個光子，先對象限空間X進行分析判斷，若進行一次四等劃分后光子在同一象限，則不進行此劃分，同時此象限劃分終止；若不在同一象限，則進行劃分，對劃分后的子象限空間按照其內(nèi)光子數(shù)迭代執(zhí)行以上兩種劃分情況，直至劃分終止。

圖6為剪枝四叉樹劃分后陰影部分光子的分隔情況，相對于圖4來說，紅框象限按照劃分準則不再進行劃分。從所對應(yīng)的樹狀結(jié)構(gòu)圖圖7可以看出，通過上述劃分方式，紅色陰影部分的枝狀結(jié)構(gòu)被剪掉，噪聲光子所處層級降低，層值變小，避免了其被誤識別為信號光子。

圖6 剪枝四叉樹分隔光子示意圖Fig.6 Photons separated by pruned quadtree

圖7 剪枝四叉樹分隔光子樹狀結(jié)構(gòu)圖Fig.7 Tree structure of photons separated by pruned quadtree

1.3.1.3 Otsu法求取光子密度閾值

由于受地表反射率、太陽高度角、大氣散射等因素的影響，沿軌方向上不同區(qū)域內(nèi)光子的噪聲密度會有所差異［11］。為適應(yīng)這種沿軌方向上信號光子和噪聲光子密度對比的不同，在沿軌方向劃分100 m等距離窗口，分窗口利用Otsu法自適應(yīng)求取各窗口光子的密度閾值，依據(jù)密度閾值進行噪聲光子和信號光子的分類。

設(shè)一個窗口內(nèi)所有光子層值取值范圍為[0，k]，假定密度閾值為t，將窗口內(nèi)的光子分為層值為[0，t-1]和[t，k]兩類，類間方差為

式中，

式中，pi為窗口內(nèi)層值為i的光子所占比，ω1、ω2分別為層值小于t和大于等于t的光子所占比，μ1、μ2分別為層值小于t和大于等于t的光子層值均值，μ為窗口內(nèi)所有光子的層值均值。

當部分噪聲光子被錯分為信號光子或者部分信號光子被錯分為噪聲光子時，就會導致類間方差變小，因此類間方差最大時意味著兩類光子錯分的概率最?。?2-23］，所以窗口內(nèi)光子的最佳密度閾值t*為

即t從(0，k)中取值，t*為使類間方差σ2(t)最大時所對應(yīng)的t值。此時窗口內(nèi)層值小于t*的光子分類為噪聲光子，其他為信號光子，依據(jù)分類結(jié)果將噪聲光子去除完成本文的一級去噪。

1.3.2 二級去噪

經(jīng)過剪枝四叉樹去噪以后，絕大部分噪點得到了消除。雖然噪聲光子在整體上比信號光子分布稀疏，但依然存在孤立噪聲光子簇，此時利用剪枝四叉樹層值表征密度時，此類局部噪聲光子由于分布密集，因此其層值偏大，被誤識別為信號光子。對于這種少量的離群噪聲光子，采用箱線圖的方法進行去除。由于箱線圖是通過光子高程四分位數(shù)來設(shè)置上下閾值進行離群點的判定，因此為顧及高程變化對箱線圖去噪的影響，同樣采用分窗口的形式進行去噪，窗口寬度在地形復雜區(qū)域和平坦區(qū)域分別取100 m和50 m。其中箱線圖的上下限為［24-25］

式中，Q1、Q3分別為窗口內(nèi)所有光子高程的下四分位數(shù)和上四分位數(shù)。如圖8為箱線圖去噪的示意圖，藍點為高程值大于Upperlimit或小于Lowerlimit的光子，其被識別為離群噪點。

圖8 箱線圖離群點檢測Fig.8 Outlier detection by box-plot

1.3.3 去噪效果驗證

為有效地驗證去噪效果，采取定性和定量相結(jié)合的方式與機載Lidar生成的DTM和DSM進行對比。定性驗證采用將四叉樹方法、本文所提方法去噪后光子擬合的地表和冠層頂部曲線，與相應(yīng)位置處DTM和DSM剖面的高程曲線進行比對的方式。在植被區(qū)地表光子和冠頂光子分別位于信號光子點云底部和頂部，由此采用移動窗口的方式，將窗口內(nèi)高程值最小和最大的信號光子分別作為地表和冠層的種子光子，然后利用三次樣條插值函數(shù)擬合得到地表和冠層曲線［9］，對比的真實地表曲線為DTM剖面高程曲線，真實冠頂曲線為DSM剖面高程曲線。定量驗證為計算種子光子高程和驗證高程數(shù)據(jù)之間均方根誤差RMSE和決定系數(shù)R2，其中RMSE和R2表示為

式中，N為計算的光子點數(shù)，hi為ICESat-2光子高程值，h?i為光子所處位置對應(yīng)的DTM或DSM高程值。

2 試驗結(jié)果與分析

本文在北達科他州研究區(qū)域1和加利福尼亞州研究區(qū)域2選擇兩條ICESat-2條帶開展去噪實驗。圖9和圖10分別為兩區(qū)域光子在沿軌方向和高程的二維分布圖，中間密集部分為信號光子，大量的噪聲光子較為稀疏地分布在信號光子帶的兩側(cè)，光子數(shù)據(jù)的信噪比皆比較低，處于強噪聲背景下。

圖9 研究區(qū)域1原始點云Fig.9 Initial point cloud in study area 1

圖10 研究區(qū)域2原始點云Fig.10 Initial point cloud in study area 2

2.1 研究區(qū)域1去噪結(jié)果

在研究區(qū)域1開展去噪實驗，圖11和圖12分別為利用四叉樹方法和剪枝四叉樹方法的去噪結(jié)果，其中四叉樹方法和剪枝四叉樹方法采用了相同的窗口劃分策略和光子密度閾值求取方法，去噪結(jié)果中藍色光子為去除的噪聲光子，紅色光子為去噪后得到的信號光子。從中，可以明顯看出，在強噪聲背景下，剪枝四叉樹方法相比四叉樹方法誤識別的噪聲光子數(shù)量得到減少，去噪效果更好。將上述兩種方法經(jīng)過二次去噪得到如圖13和14所示結(jié)果圖，可以看到四叉樹方法去噪結(jié)果在經(jīng)過二次去噪后仍存在一定量的噪聲光子，而本文所提剪枝四叉樹方法和二級去噪則基本去除了所有噪聲光子。這是由于二次去噪所采用的箱線圖法的去噪性能會隨著離群噪點的增多而降低，而剪枝四叉樹方法減少了誤識別噪聲光子的數(shù)量，所以使剪枝四叉樹結(jié)合二級去噪在去噪效果上要優(yōu)于四叉樹方法結(jié)合二級去噪。

圖11 四叉樹方法在研究區(qū)域1去噪結(jié)果Fig.11 Denoising result of quadtree method in study area1

圖12 剪枝四叉樹方法在研究區(qū)域1去噪結(jié)果Fig.12 Denoising result of pruned quadtree method in study area1

圖13 四叉樹方法和二級去噪在研究區(qū)域1去噪結(jié)果Fig.13 Denoising result of quadtree method and second-level denoising in study area 1

圖14 剪枝四叉樹方法和二級去噪在研究區(qū)域1去噪結(jié)果Fig.14 Denoising result of pruned quadtree method and second-level denoising in study area 1

為進一步驗證本文所提方法相較于四叉樹方法具有更優(yōu)的去噪性能以及去噪結(jié)果的準確程度，本文以機載激光雷達生成的高程數(shù)據(jù)產(chǎn)品為標準進行對比試驗。由于此區(qū)域植被覆蓋稀疏，因此在研究區(qū)域1只使用DTM作為對比數(shù)據(jù)。分別用四叉樹方法提取的信號光子和本文所提方法提取的信號光子得到地表種子光子，將擬合的地表曲線與DTM剖面高程曲線進行對比，結(jié)果如圖15所示。可以看到，利用四叉樹方法提取的信號光子所擬合的地表曲線與DTM剖面的真實地表曲線之間存在著較大偏差，無法正確反映地表走勢，而本文所提方法提取的信號光子擬合的地表曲線則與DTM剖面的地表曲線基本一致。進一步地，定量計算兩種方法所得地表種子光子高程與DTM高程之間的均方根誤差與決定系數(shù)，其結(jié)果如表1所示。可以看出，本文方法的RMSE為0.91 m，優(yōu)于四叉樹方法所對應(yīng)的10.51 m，R2為0.997，非常接近理想值1，計算光子的高程與分辨率為1 m的機載激光雷達高程數(shù)據(jù)產(chǎn)品高程基本一致，較四叉樹方法對應(yīng)的0.713，有顯著提高。

圖15 研究區(qū)域1地表曲線結(jié)果Fig.15 Ground curves in study area 1

表1 研究區(qū)域1地表精度評價Table 1 Ground accuracy evaluation in study area 1

2.2 研究區(qū)域2去噪

圖16和圖17分別為利用四叉樹方法和剪枝四叉樹方法在研究區(qū)域2的去噪結(jié)果，可以看出，四叉樹方法誤識別的噪聲光子較多，而剪枝四叉樹方法明顯減少了誤識別噪聲光子的數(shù)量。圖18為四叉樹方法和二級去噪的結(jié)果，仍然殘留了一定量的噪聲光子，圖19所示的剪枝四叉樹方法和二級去噪則基本去除了所有的噪聲光子。

圖16 四叉樹方法在研究區(qū)域2去噪結(jié)果Fig.16 Denoising result of quadtree method in study area 2

圖17 剪枝四叉樹方法在研究區(qū)域2去噪結(jié)果Fig.17 Denoising result of pruned quadtree method in study area 2

圖18 四叉樹方法和二級去噪在研究區(qū)域2去噪結(jié)果Fig.18 Denoising result of quadtree method and second-level denoising in study area 2

圖19 剪枝四叉樹方法和二級去噪在研究區(qū)域2去噪結(jié)果Fig.19 Denoising result of pruned quadtree method and second-level denoising in study area 2

同樣為進一步驗證本文所提方法相較于四叉樹方法具有更好的去噪性能以及去噪結(jié)果的準確程度，對于研究區(qū)域2，從四叉樹方法提取的信號光子和本文所提方法提取的信號光子得到冠頂和地表種子點，將擬合的冠頂和地表曲線分別與DSM和DTM剖面高程曲線進行對比，如圖20和21所示。從圖中可以看出，利用四叉樹方法提取的信號光子擬合的冠頂和地表曲線與DSM和DTM剖面的冠頂和地表曲線都存在較大的偏差，而利用本文所提方法提取的信號光子擬合的曲線則與DSM和DTM剖面曲線基本一致。從表2和表3可以看出，研究區(qū)域2本文方法所得信號光子中提取的冠層和地表種子光子高程與DSM和DTM高程之間的RMSE分別為3.56 m和2.47 m，明顯優(yōu)于四叉樹方法對應(yīng)的90.17 m和90.92 m，R2為0.998和0.999，接近理想值1，計算光子的高程與分辨率為1 m的機載激光雷達高程數(shù)據(jù)產(chǎn)品高程基本一致，與四叉樹對應(yīng)的0.4和0.156相比，有顯著提高。

圖20 研究區(qū)域2冠頂曲線結(jié)果Fig.20 Canopy top curves in study area 2

圖21 研究區(qū)域2地表曲線結(jié)果Fig.21 Ground curves in study area 2

表2 研究區(qū)域2冠頂精度評價Table 2 Canopy top accuracy evaluation in study area 2

表3 研究區(qū)域2地表精度評價Table 3 Ground accuracy evaluation in study area 2

3 結(jié)論

本文提出了一種無輸入?yún)?shù)強噪聲背景下ICESat-2點云去噪方法，以剪枝四叉樹層值表征光子密度，自適應(yīng)求取信號光子密度閾值，進行一級去噪；然后通過箱線圖進一步去除少量離群噪聲光子，進行二級去噪。本文所提方法避免了強噪聲背景下局部稀疏、相距較近的噪聲光子被四叉樹分隔層次較多而被誤表征為密度過大，自適應(yīng)求取信號光子的密度閾值，并對相聚密集噪聲光子進行了有效去除，與四叉樹方法相比，具有更優(yōu)的去噪性能。本文所提方法去噪后的信號光子所擬合的冠頂和地表曲線與分辨率為1 m的機載激光雷達高程數(shù)據(jù)產(chǎn)品高程剖面曲線基本一致。

當然，限于所掌握的公開數(shù)據(jù)有限等原因，本文只選取了較為典型的區(qū)域進行了比對實驗，更多區(qū)域的實驗驗證還有待以后進一步進行。同時，本文去噪方法中采用了劃分窗口的方式，可能會存在不同窗口去噪結(jié)果相鄰邊界不一致的現(xiàn)象，后續(xù)研究中將結(jié)合窗口拓展等方式來降低此類邊緣效應(yīng)。