王濱輝：全波形激光雷達的波形優化分解演算法

《測繪學報》

構建與學術的橋樑 拉近與權威的距離

全波形激光雷達的波形優化分解演算法

王濱輝1,2, 宋沙磊2, 龔威1, 陳振威2, 林鑫2, 程學武2, 李發泉2, 史碩1

1. 武漢大學測繪遙感信息工程國家重點實驗室, 湖北 武漢 430079;

2. 中國科學院武漢物理與數學研究所波譜與原子分子物理國家重點實驗室, 湖北 武漢 430071

收稿日期：2017-01-24；修回日期：2017-09-08

摘要：隨著數據存儲能力和處理速度的提高，三維激光掃描系統逐漸具備全波形採集和分析技術。為了從全波形數據中獲得脈衝時間、幅度、脈寬以及多回波分布等綜合信息，波形分解成為了全波形激光雷達數據處理的關鍵技術之一。針對LM演算法在一定程度上依賴初值，而傳統激光雷達數據處理容易遺漏部分重疊的返回波，本文提出了一種改進回波分量初值設定的演算法來獲取回波脈衝的位置、寬度和強度。針對一套自主研發的全波形記錄激光雷達演示系統進行了波形分解試驗，定性和定量分析結果驗證了該方法的有效性、可靠性和準確性。

關鍵詞：全波形LiDAR波形分解高斯函數LM演算法

Optimization Decomposition Method of Full-waveform LiDAR

WANG Binhui1,2, SONG Shalei2, GONG Wei1, CHEN Zhenwei2, LIN Xin2, CHENG Xuewu2, LI Faquan2, SHI Shuo1

Abstract: With the improvement of data storage capacity and data-processing capabilities, full waveform LiDAR develops rapidly and then from waveform data abundant information about the physical characteristics of the targets can be effectively retrieved through data processing. Waveform decomposition is therefore the key to full waveform LiDAR data processing. However, the number and initial parameters of the echo components are difficult to set in waveform decomposition. Conventional decomposition methods detect echo components by peak points or using the threshold method, which may ignore some overlapping components and show low accuracy. To this end, in this paper a novel method based on LM algorithm that takes into account both peak points and inflection points is adopted and it can extract the location, amplitude and FHWM of the echo components, proving it is a reliable and high accurate decomposition algorithm. To further demonstrate the advantages of the suggested method, waveform data measured by a full waveform LiDAR demonstration system and generated from simulation were both decomposed using the method. The results show that the suggested algorithm is efficient, promising and can effectively decompose adjacent echo components, then will improve the accuracy in the next phase of data processing.

Key words:full-waveform LiDARwaveform decompositionGaussian functionLM algorithm

機載激光雷達(light detection and ranging, LiDAR)集全球定位系統、慣性導航系統和激光掃描系統於一體，是一種快速獲取地表三維信息的主動式探測技術[1]。早期的機載激光雷達能夠記錄有限個離散的返回信號，通過發射脈衝和回波脈衝之間的時間間隔獲取距離。隨著數據存儲能力和處理速度的提高，新一代的三維激光掃描系統逐漸具備全波形採集和分析技術。近十多年來，全波形激光雷達系統發展迅速，被廣泛應用於地基、機載、車載和星載等平台。1999年，NASA設計了機載激光植被成像感測器[2](laser vegetation imaging sensor, LVIS)；2003年，NASA發射了搭載激光測高系統[2](geoscience laser altimeter system, GLAS)的ICESAT(ice, cloud and land elevation satellite)衛星；國內外對機載小光斑全波形激光雷達數據處理的相關研究起始較晚。2004年，RIEGL公司首次在商業LiDAR系統中提供了具備全波形數字化能力的小光斑激光雷達掃描儀[3]。目前，全球主要的商業全波形激光雷達系統主要有瑞士的TopEye MarkⅡ、奧地利的Rigel LMS-Q560、LMS-Q680和加拿大的Optech Inc ALTM3100[4]。全波形激光雷達能夠以很小的採樣間隔記錄散射體的整個後向散射回波波形，其中所蘊含的豐富信息量很快引起了本研究領域學者的極大興趣，開啟了全波形激光雷達數據處理與應用研究時代。

相對於離散激光雷達系統，全波形激光雷達系統能夠提供更多的目標信息[5]，單次激光脈衝可獲取一個激光腳印內複雜地物目標的回波形態，可充分分析脈衝時間、幅度、脈寬及多回波分布等綜合信息[6]。波形分解演算法對地物目標的垂直結構，尤其是植被冠層的精細結構分析[7]有很大的幫助。因此，全波形數據的處理在大氣監測、三維城市建模、地形測繪與可視化、機器人導航和電力線檢測等民用領域和戰場偵察、精確制導等軍事領域都得到了廣泛應用。目前，已有很多演算法應用到波形處理研究中，通常將其分為3類[3]：閾值法、反卷積法和波形分解法。閾值法對脈衝個數較少的波形處理效果好，回波脈衝的探測能力有限。反卷積法[8]通過將已知的發射脈衝和回波信號作反卷積處理來獲取地物的後向散射回波信息。與閾值法相比，這種方法雖然提高了回波檢測能力，卻仍然沒有分解出波形參數[9]。波形分解法認為回波信號是激光路徑內各個地物對激光脈衝響應的總和。文獻[10]提出了波形數據可以用一系列高斯函數擬合。文獻[11]考慮到回波分量波形的多樣性，進一步提出了廣義高斯模型，並使用Levenberg Marquardt(LM)演算法[12]進行參數優化。文獻[13]提出了用於探測和定位全波形數據中回波分量的均方差方法。文獻[14]將期望最大化演算法(expectation maximation, EM)用于波形擬合。文獻[15]考慮波形分解的物理意義和數學意義，提出用改進的EM演算法優化模型參數。文獻[16]提出逐級遞進分解策略解決高斯函數個數的確定問題。文獻[17]針對微弱回波分解能力不足的缺點，提出一種迭代的波形分解方法。

傳統LM演算法[18-19]在一定程度上依賴初值，而激光雷達數據處理一般採用固定數的波形分解方法[20]，容易遺漏部分重疊的返回波。針對波形分解的研究現狀，本文提出一種全波形激光雷達的優化分解演算法(optimization decomposition method, ODM)，該演算法通過改進回波分量初值設定來獲取回波脈衝的位置、寬度和強度[21-22]，在依次輸入初始回波分量的同時，利用LM演算法優化模型參數，剔除小於誤差閾值的分量。針對一套全波形激光雷達演示系統獲得的大量全波形數據進行了波形分解試驗，定量分析結果驗證了該方法對於分解相近回波分量的有效性、可靠性和準確性。

1 全波形記錄原理

全波形激光雷達的波形，是指激光發射脈衝或接收脈衝的能量隨著時間變化的函數。全波形激光雷達系統可對同目標接觸和相互作用後反射回來的後向散射回波脈衝進行小時間間隔的強度採樣，對強度進行數字量化，並記錄下後向散射脈衝的強度值。後向散射回波波形是激光脈衝所照射的激光足印中所有目標與激光脈衝接觸、相互作用並後向散射後的能量與系統雜訊的疊加效果[4]。全波形激光雷達對發射脈衝和後向散射回波脈衝都進行小間隔採樣，獲取的原始數據如圖 1所示，幾乎能記錄整個後向散射回波波形。

與傳統激光雷達系統記錄點雲數據不同，全波形激光雷達系統能夠以波形的形式記錄一定高程範圍內不同高程點上的後向散射能量[23-24]。理論上回波波形可以看作是若干個高斯函數的疊加，如圖 1(b)所示，高斯函數的振幅、峰值位置、波形寬度、峰值與峰值之間的間距等都是波形分析中的重要參數[25-27]。

設理想波形模型為高斯混合模型，表達為k個高斯函數分量的疊加

(1)

式中，Ai、μi、Fi分別表示第i個高斯分量的脈衝振幅、中心位置和半高寬；noise是原始波形的背景雜訊。波形擬合實質是通過n個採樣點，逐步迭代計算得到3×k個高斯參數和一個背景雜訊參數的最優解，使得擬合波形與原始波形相差儘可能小。

2 全波形優化分解演算法

針對易遺漏的重疊的回波分量，本文提出了改進的預估回波分量個數和回波分量初值的演算法。圖 2所示為該演算法總體流程，其中，實線框為主要分解步驟，虛線框為相應步驟的補充說明。為了詳細描述波形分解演算法的處理過程，選取圖 1(b)作為實例。這段波形基於Matlab生成，含有5個高斯分量，544個採樣點。高斯分量的振幅、中心位置和半高寬參數見表 1。

表 1實例波形的模擬參數Tab. 1Waveform simulation parameters

表選項

2.1 全波形數據預處理

通常，為保證波形記錄沒有遺漏，所記錄的波形範圍會超過要求的範圍。於是，回波數據的開頭部分或者末尾部分會用來估算雜訊水平。開頭部分或結尾部分的範圍依試驗數據的特點而定，盡量取不存在回波分量的背景雜訊。μn和σn分別表示雜訊均值和標準差。探測目標的距離範圍不定可能造成開頭部分或末尾部分有回波分量的存在，於是雜訊水平可能被高估。考慮到開頭和結尾同時存在回波分量的可能性很低，選擇其中一個部分，如圖 3所示，來進行雜訊水平的估算。首先，μf、μr、σf、σr分別表示前後部分的雜訊均值和前後部分的標準差。如果μfμr，則μn=μf，σn=σf；否則，μn=μr，σn=σr，雜訊閾值thn計算公式為

(2)

波形在採集的過程中由於多方面的原因會產生背景雜訊，這些雜訊一般呈現小振幅抖動。原始波形中的雜訊可以被看作是高斯白雜訊，回波分量採用高斯函數混合模型擬合，高斯平滑濾波器可以實現原始波形的最優匹配濾波[28]。濾波器的帶寬設為發射波的半高寬。

2.2 全波形參數初始化估計

原始波形數據經過高斯平滑濾波器預處理之後，大部分隨機雜訊被消除，處理後的波形更接近真實回波信號，有利於下一步的波形分解。為了不遺漏回波分量，本文除了採用一般的波形分解方法即利用局部峰值點確定回波分量之外，還進行了二次峰值點探測，確定重疊的回波分量。

2.2.1 全波形分量初判斷

計算每個採樣點的一階導數，即一階差分，卷積核為[-1, 0, 1]，將零交叉值(由正變負)的點視為波峰位置；計算每個採樣點的二階導數，即二階差分，卷積核為[1, -2, 1]，將零交叉值的點視為拐點位置，剔除小於雜訊閾值thn的峰值點和拐點。

對於每一個峰值點，估算它所在回波分量的初值。對於峰值點Pi，若其左右單調區間內均無拐點，則認為Pi不屬於任何回波分量。否則，分別計算其左右單調區間內拐點位置的均值，記為μiL、μiR。分別計算μiL、μiR與Pi的距離，記為diL、diR。若diL、diR均小於發射波半高寬dFWHM/2，則認為Pi不屬於任何回波分量。否則，Pi可能屬於某一回波分量，估算其初值。將Pi的振幅作為脈衝振幅初始值Ai，Pi的位置作為中心位置初始值μi，半高寬初始值Fi計算公式為

(3)

式中，半高寬閾值為dFWHM。當左右區間內拐點平均位置到峰值點距離diR、diL均不小於dFWHM/2時，取diR、diL較小值的2倍作為Fi的初始值。當diR、diL中的一個不小於dFWHM/2，另一個小於dFWHM/2時，取diR、diL較大值的2倍作為Fi的初始值。

初判斷出來的回波分量的峰值點如圖 4所示，圖中找出了8個可能存在的回波分量所對應的峰值點。

2.2.2 全波形分量二次判斷

根據峰值點檢測回波分量，大部分距離較遠或重疊度不高的回波分量被探測了出來。然而，這種方法很容易遺漏重疊的回波分量[29]。針對這個問題，進行二次探測回波分量。

對於初判斷中的每個回波分量，搜索其峰值點左右區間內的拐點，若存在不屬於區間[μ-Fi/2,μ+Fi/2]的連續拐點，且這組拐點間的最大距離大於dFWHM/2，則說明附近可能存在其他的回波分量。將這組連續拐點中的振幅較大值作為新的回波分量的脈衝振幅初始值Aj，將具有較大振幅的拐點的位置作為新的回波分量的中心位置初始值μj，將這組拐點間的最大距離的2倍作為新的回波分量的半高寬初始值Fj。如圖 5所示，經過二次判斷，又找出了1個可能存在的回波分量所對應的峰值點。

2.3 參數最優化估計

將初始化估計的回波分量按能量大小降序排列，回波能量用振幅與半高寬的乘積(Ai×Fi)表示。採用LM演算法，按順序向波形分解函數模型中輸入初始波形參數，迭代計算最優解。每次添加一組波形參數，就更新一次波形分解模型參數，剔除其中幅值小於雜訊閾值或者半高寬小於發射波半高寬dFWHM的回波分量，併合並相對距離小於dFWHM/2的兩個初始回波分量。直到原始波形與擬合波形之間的殘差小於雜訊的標準差σn或者波形分解參數達到上限Nmax，停止輸入初始波形參數。

3 波形分解實測及模擬驗證3.1 實測驗證與分析

本文使用的數據來自一套全波形激光雷達演示系統，該系統各項參數見表 2。試驗時，激光垂直照射到兩塊白板上，兩塊白板平分一塊光斑。兩塊白板相距10、20、30、40、50、60、70、80、90、100 cm。每種間隔選取1000個回波信號進行波形分解，提取各項波形參數。每段回波信號的採樣點個數為30~40個。對於實測回波信號，取開頭10個採樣點和結尾10個採樣點用作雜訊評估。

表 2激光雷達演示系統參數Tab. 2Specifications of the LiDAR demonstration system

表選項

由於實際測量時波形參數的真實值未知，於是通過決定係數R2來評價波形擬合的精度

(4)

式中，si表示原始波形中第i個採樣點的強度值；yi表示波形分解模型中第i個點的強度值。決定係數R2是衡量兩個隨機變數之間線性相關程度的指標，決定係數越接近1，說明擬合數據與原始數據相關性越強，波形擬合的精度越高。雖然決定係數R2可以評價擬合波形整體與原始波形的近似程度，但不能說明分解相近回波分量的能力。於是計算分解得到的波形分量總數與真實波形分量總數的比值Rt。比值Rt越接近1，說明分解的波形分量個數越接近真實的目標個數，分解相近回波分量的能力越高。

本文對相對距離分別為10、20、30、40、50 cm的5組原始數據採用ODM方法和傳統LM演算法進行了對比分析，如圖 6所示。當兩回波分量的相對距離是10 cm時，如圖 6(a)所示，採用兩種方法的分解效果相近，都從原始數據中分解出一個回波分量，並且有較高的擬合決定係數R2值(大於0.98)和較小的擬合殘差RMSE值(2左右)。當兩回波分量的相對距離是20、30、40 cm時，如圖 6(b)-(d)所示，雖然採用兩種方法進行波形分解都有較高的擬合決定係數R2值和較小的擬合殘差RMSE值，即擬合的波形都能很好地近似原始波形，但是使用ODM方法可以分解出兩個回波分量，而傳統的LM演算法只能分解出一個回波分量，這說明ODM方法分解相近回波分量的效果更好，能夠更好地從原始數據中提取出更多的目標信息。當兩回波分量的相對距離是50 cm時，如圖 6(e)所示，採用兩種方法的分解效果相近，都從原始數據中分解出兩個回波分量，並且有較高的擬合決定係數R2值(大於0.99)和較小的擬合殘差RMSE值(小於2)，這說明對於相對距離超出某個範圍的兩個回波分量，採用兩種方法得到的分解效果相近。

對選取的10 000個回波信號的分解結果進行統計分析，分解得到的波形分量總數與真實波形分量總數的比值Rt以及決定係數R2的平均值和標準差見表 3。由表 3可知，ODM演算法成功分解出85.75%的回波分量，遠高於採用傳統LM演算法分解出的波形分量數60.08%，這說明ODM分解相近回波分量的效果更好，能夠更好地從原始數據中提取出更多的目標信息。採用ODM演算法分解結果的擬合決定係數均值和標準差分別是0.962 2和0.107 6，採用傳統LM演算法分解結果的擬合決定係數均值和標準差分別是0.870 3和0.301 5，這說明採用ODM演算法擬合的波形更接近原始波形，擬合精度更高。

表 3波形分解結果統計(實測數據)Tab. 3Statistics on waveform decomposition results (measured data)

表選項

3.2 模擬測試與分析

為進一步驗證回波分量的相對距離對波形分解可行性的影響，利用Matlab模擬了只含兩個回波分量的波形數據。模擬數據的波形參數(振幅、中心位置、半高寬)是一系列隨機整數，其範圍見表 4。模擬回波信號開頭50個採樣點和末尾50個採樣點被用作評估雜訊水平。模擬輸入的各項波形參數作為真實值，可與分解結果對比來定量評估波形參數的精確性。以振幅參數為例，AT表示振幅的真實值，AC表示振幅的優化結果，則振幅的偏差EA

(5)

表 4模擬數據波形參數列Tab. 4Simulation data parameters

表選項

分別計算EA的均值μeA和標準差σeA，這兩個誤差評估參數是評價振幅的可靠性和準確性的重要指標。同樣，中心位置和半高寬兩項波形參數的誤差評估參數記為μeL、σeL、μeF、σeF。對5000個模擬信號的分解結果進行統計分析，其分解得到的波形分量總數與真實波形分量總數的比值Rt及各項誤差評估參數的平均值和標準差如表 5所示。由表 5可知，雖然採用兩種方法分解得到各項波形參數的誤差水平相近，但是ODM演算法成功分解出83.75%的回波分量，遠高於採用傳統LM演算法分解出的波形分量數55.73%，這說明ODM分解相近回波分量的效果更好，能夠更好地從原始數據中提取出更多的目標信息。

表 5波形分解結果統計(模擬數據)Tab. 5Statistics on waveform decomposition results(simulation data)

表選項

4 結論

相對於離散激光雷達系統，全波形激光雷達系統能夠提供更多的目標信息。發射單次激光脈衝可獲取一個激光腳印內複雜地物目標的回波形態，可充分分析包括脈衝時間、幅度、脈寬以及多回波分布等全波形綜合信息。傳統LM演算法在一定程度上依賴初值，而激光雷達數據處理一般採用固定數的波形分解方法，容易遺漏部分重疊的返回波。針對波形分解的研究現狀，本文提出ODM方法，該演算法通過改進回波分量初值設定來獲取回波脈衝的位置、寬度和強度，在依次輸入初始回波分量的同時，利用LM演算法優化模型參數，剔除小於誤差閾值的分量。利用ODM方法對實測試驗數據和模擬試驗數據進行了波形分解試驗，定性和定量地驗證了該方法的有效性、可靠性和準確性。ODM方法不僅能很好地擬合原始波形，還能有效地分解相近的回波分量，能夠更好地從原始數據中提取出更多的信息。對於一些非對稱或拖尾的波形，該方法還不能非常準確地擬合。在未來的工作中，將探索如何提高波形參數的精度，提高分解相近回波分量的能力，使演算法更普遍地適用於各種形狀的波形。例如，可以考慮用其他隨機模型代替高斯模型來擬合波形分量，從而提高波形擬合的精度。

【引文格式】王濱輝，宋沙磊，龔威，等。全波形激光雷達的波形優化分解演算法[J]. 測繪學報，2017，46(11)：1859-1867. DOI: 10.11947/j.AGCS.2017.20170045

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！