機載小光斑全波形LiDAR數據處理及應用

周靜平，張愛武，王書民

（首都師范大學三維信息獲取與應用教育部重點實驗室，北京 100048）

一、引 言

機載激光雷達（airborne laser scanning，ALS）也指雷達或激光雷達，是一種通過激光器向地表發射紅外線短波脈沖，并由光電接收器接收背向散射的回波信號的主動式光學遙感器［1］。作為一種先進的新型對地觀測系統，它集激光、全球定位系統和慣性導航系統（inertial navigation system，INS）3項技術于一身，并用于快速獲取地面目標的三維空間信息。

自20世紀70年代美國、加拿大率先研制航空激光雷達測量系統以來，此項技術發展特別迅速［2］。20世紀80年代開始出現一些試驗系統，到90年代中后期逐漸成熟。尤其是20世紀90年代中期商業ALS產品和服務制造商的出現，使得激光雷達系統和存儲技術得到了前所未有的迅猛發展。2004年，奧地利Riegl公司推出了世界上第一臺商業機載小光斑全波形激光雷達數字測量系統。隨之，Optech ALTM 3100E、TopEye Mark II、Leica ALSII等機載激光雷達系統也均可進行全波形數字化測量，并在地形測繪、林業、環境監測、三維城市建模等許多領域得到了廣泛的應用［3］。

二、全波形LiDAR理論

1.激光傳感器

激光雷達的工作原理與傳統雷達遙感大體相同。激光測距系統是其核心組成部分，以激光波束掃描的方式測量激光器到目標方向激光照射點的距離，即通過測量地面點的回波脈沖與激光的發射脈沖之間的時間延遲得到激光器到地面點之間的距離［4］。測距原理公式為

式中，R為激光器到目標物的距離；c為光速；t為激光脈沖從激光器到目標物傳輸的往返時間［5］。機載激光雷達系統的R和t是可以相互轉換的，具體使用哪個可根據實際需求進行選擇。激光雷達的測量方式是十分精確的，其測量精度常?？蛇_到毫米級。再結合GPS記錄的坐標信息和INS記錄的姿態角信息，通過聯合解算便可得到激光點的大地坐標。

2.全波形激光雷達

全波形LiDAR即全數字波形激光雷達，采用數字化方式，不僅記錄若干次離散回波信號，而且將發射信號和回波信號均以很小的采樣間隔進行采樣并記錄［6-8］。激光雷達記錄的回波波形是對光斑內各個點的反射信號按時間先后順序的記錄，回波波形代表了激光后向散射的能量，也就是說，回波波形可以看做是回波強度信息在接收時間軸上的一個函數。

較之傳統的LiDAR，全波形LiDAR具有3個方面的優點［7］:① 通過對接收的回波信號進行處理，可以提取出所有的脈沖回波信息，也就是說同一激光束中，全波形LiDAR系統比傳統LiDAR系統能提供更詳細的點云信息（特別是在林木覆蓋區域），且密度大，層次感強；②通過對回波進行波形分解、擬合，能得到更加精確的電云坐標；③ 對接收到的波形信號進行處理、建模，通過模型可得到更加豐富的地物特征。全波形LiDAR對波形信號進行了精確細致的存儲和刻畫，尤其是對于樹木等特殊地物，它能提供更多的目標信息。

3.雷達方程

結合傳統雷達方程，激光雷達方程可表示如下［1］

式中，Pr是接收的能量；Pt是發射脈沖能量；Dr為接收器孔徑大??；R是目標與飛行器的距離；βt為發散角；ηsys表示系統衰減系數；ηatm表示大氣的衰減系數；σ被稱為背向散射截面（back scatter cross section）。很顯然，式（2）是在理想狀況下的一個簡單散射體的散射能量方程，對于一個光斑內存在多個目標物的復雜散射體，其反射的能量組成將更為復雜，仍需進一步研究。

波形通常被模型化為一系列高斯脈沖。根據下面公式進行推算，最終得到接收的能量Pr被表示為［1］

式中，ti為往返的時間；為目標i的振幅；Sp，i為回波脈沖的標準差。

式中，Ss是發射脈沖的標準差；為發射脈沖的振幅；Si是目標i回波的標準差；σi為目標i總的后向散射截面，是目標地物不同散射截面的積分。目標的散射截面與地表的粗糙度有關，可根據發射脈沖信號與不同目標的后向散射截面的積分得到一個接收的波形。

對全波形數據進行處理時，通常還要進行輻射校正。本文采用的是Riegl LMS-Q560數據，經過分析比較后決定采用Wolfgang Wagner所提出的利用柏油路面（反射率為0.2）作為標準參考的方法，具體校正公式為［9］

三、波形分解方法和步驟

1.分解方法

波形分析技術最先由美國宇航局（NASA）用于其機載大光斑LiDAR系統上，激光雷達波形數據分析也是以美國為首進行發展的。目前，對機載激光全波形數據進行分解，國外研究得較早也較為深入，國內研究最近幾年也逐漸開始。

目前波形數據的處理方法主要有3種:①閾值法 （average square difference function［10］，ASDF）；②反卷積法；③波形分解法。

波形分解的方法是將回波的波形看作是一個光斑內不同高度的目標物對激光脈沖反射后綜合作用［11］的結果，其將不同高度處的目標物回波從接收的波形信號中提取出來，從而達到對目標進行定位和識別的目的。波形分解的核心內容包括噪聲提取和目標物反射回波波形的準確定義。目前波形分解的研究都是基于對回波理想情況下的模擬，多采用高斯函數或者正態分布函數近似刻畫單個目標的回波波形。其中，采用高斯函數的形式是目前采用較多也是較為準確的一種方法。Wagner［12］等最先采用了高斯分解的方法對第一個能夠記錄完整波形小光斑激光雷達Riegl LMS-Q560的數據進行了處理，結果證明98%以上的回波數據可以用高斯函數很好地擬合，其精度也能夠達到規定的要求。目前，有幾種方法可以實現高斯分解，如Wagner等開發的基于非線性最小二乘的分層高斯函數波形擬合法；Jutzi等實現的高斯-牛頓法；Persson等的 EM（Expectation-Maximization）［13］算法實現；武漢大學的馬洪超、李奇［14］采用改進的EM算法的波形分析技術；De Boor和Cox分別提出的遞推定義的B樣條法等。波形分解是目前機載激光雷達［15-16］較為常用的波形數據處理算法，均可獲得較為理想的結果，并可獲取更多的額外參數。

本文采用基于高斯數學模型的波形分解法，對研究區波形數據進行波形分解和擬合，采用L-M算法［1］并利用最小二乘法原理不斷迭代優化，最終得到最優的參數，即每個單個波形的振幅、距離、波寬。

2.分解步驟

1）數據預處理（波形數據讀取、噪聲去除、波形平滑）。

2）參數初始化（峰值點檢測與初始參數估計）。

3）參數提?。ㄒ罁M合函數和分解算法）。

四、試驗成果和分析

本文的研究區位于河北省城區，該區域內地物類型豐富，有人工建筑物、道路、高大樹木、低矮灌木、草地等。經過對研究區機載小光斑全波形Li-DAR數據進行波形分解，得到了5 481 859條波形數據。為了直觀地顯示及后續處理的方便，本文在Matlab平臺上，制作出了波形數據的顯示界面，如圖1所示。

圖1 波形顯示界面

利用此界面，可對本研究區域所有的波形數據逐條地進行直觀的顯示，也為后期的數據處理提供了很多便利。研究區具有豐富的地物類型，所得的數據回波信息也很豐富，回波次數最多可達7次，詳細情況見表1。

表1 回波數目分布情況

在Matlab平臺上，對研究區全波形數據進行了顯示，并且利用激光點云數據處理專業軟件Terrasolid對點云數據進行了大致分類（見表2）。由于研究區域較大，為了方便顯示，只截取了其中一部分，效果如圖2所示。

表2 分類情況

從表1和圖1可以明顯看出，研究區的回波次數最多為7次，最少為1次。大部分為1次回波占87.814%；2次以上回波僅占12.186%。這是因為研究區位于城區，雖然地物類型較為豐富，但因建筑用地（房屋建筑、道路等）面積占據了大部分（占86.92%），這些地區較為平坦，激光對其無穿透；此研究區有也較多高大樹木、低矮灌木等（占13.08%），林木等地物具有葉間、枝間縫隙，激光可穿透。對比地物分類圖，可以明顯看出，建筑用地及道路等非植被區域和一次回波區大致重合，僅有1次回波；而有林木分布的植被區域和多次回波區大致吻合，回波為多次。整體上，也可以從表1和表2中1次回波的百分比和非植被區的百分比大致相同，多次回波的總百分比和植被區的百分比大致相等來進一步得到說明。值得注意的是，研究區還有部分回波為0，這是因為地物反射波的方向不在激光接收器的視域內，導致接收不到地物的回波信息，或者是地物的反射率太低所致。

圖2 回波次數及地物分類圖

另外，第一次回波反映的是建筑物輪廓信息和高大灌木的冠層，可用于道路信息、建筑物的輪廓信息提取，也可用于計算高大灌木的植株高度信息；其他次回波反映的是林木信息，可用于計算林木的蓄積量；最后一次回波反映的是地面，可用于輔助地形圖和數字高程模型制作。

經過波形分解，提取出了回波脈沖的距離、振幅和波寬等參數信息，結合點云分類信息，所成圖件如圖3所示。

研究區波形數據所提取的參數信息中，距離參數指的是地物與激光發射器之間的距離，間接反映了高程上的差異，傳統的分類有很多都是基于高程進行的；振幅參數與地物類型間關系最為緊密，反映也最為直接；波寬參數與植被間關系較為緊密。該研究區所提取出的波形數據，振幅最大值為0.844 3 m，最小值為0.000 137 m，大多分布在0.06 m以下（即0 m～0.01 m），如圖4所示；波寬最大值為214 nm，最小值為14 nm，大多分布在150 nm以下。

圖3 回波參數及分類圖

值得說明的是，距離、振幅、波寬信息只是波形的單一參數，若僅僅用它進行分類，效果并不十分理想。由圖3可知，植被對波寬參數比較敏感，并且隨著樹木高度的增加，其波寬也隨之增大。如高大樹木比低矮草地的波寬要寬；高大建筑物和樹木冠層對距離信息比較敏感。上述距離、振幅、波寬是波形數據的重要參數信息，可以為后期點云精分類和數據的深入應用提供可靠依據。因此，這些參數信息的提取是十分重要和必要的。

圖4 振幅分布圖

五、結論與展望

本文以高斯函數波形分解理論為依據，對研究區的機載小光斑全波形LiDAR數據進行了波形處理，提取出了各個波形的振幅、距離、波寬等參數，這些特征對后續的地物分類提供了科學依據。通過試驗，充分驗證了利用高斯模型函數來擬合機載小光斑全波形脈沖是切實可行的，并在Matlab平臺上開發出了機載小光斑全波形數據顯示界面，可方便直觀地顯示所想查看的各條波形信息。有效波形的提取是波形分析領域里的一個難點問題。此外，背向散射體截面是一個很主要的參數，它包含了不同地物的距離特征和散射特征，這個參數對點云的分類指導意義特別重大。如何對機載小光斑全波形數據進行有效的波形分解和如何充分利用分解后提取出的參數信息對地物進行分類，將是進一步研究的內容。

機載小光斑全波形數據具有能夠快速獲取地表物體的精確三維空間信息和很多的地物特征信息的優勢，未來還將繼續在地形測繪、農林業、數字城市建設等諸多領域得到越來越廣泛和深入的應用。

［1］ WOLFGANG W.Gaussian Decomposition and Calibration of a Novel Small-footprint Full-waveform Digitizing Airborne Laser Scanner［J］.ISPRS Journal of Photogrammetry ＆ Remote Sensing，2006，60:100-112.

［2］ ACKERMANN F.Airborne Laser Scanning-Present Status and Future Expectations［J］.ISPRS Journal of Photogrammetry ＆ Remote Sensing，1999，54（2-3）:64-67.

［3］ 王成，MENENTI M.機載激光雷達數據的誤差分析及校正［J］.遙感學報，2007，11（3）:390-397.

［4］ BACHMAN C G.Laser Radar Systems and Techniques［M］.London:Artech House，1979.

［5］ BALTSAVIAS E P.Airborne Laser Scanning:Basic Relations and Formulas［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1999，54（2-3）:199-214.

［6］ HOFTON M A，MINSTER J B，BLAIR J B.Decomposition of Laser Altimeter Waveforms［J］.IEEE Transactions on Geoscience and Remote Sensing，2000，38（4），1989-1996.

［7］ HUG C.A Waveform Digitizing LiDAR Terrain and Vegetation Mapping System［J］.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences（Freiburg，Germany），2004（36）:24-29.

［8］ PERSSON á，S?DERMAN U，T?PEL J，et al.Visualization and Analysis of Full-waveform Airborne Laser Scanner Data［J］.International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2005（36）:103-108.

［9］ WAGNER W，RONCAT A，MELZER T，et al.Waveform Analysis Techniques in Airborne Laser Scanning.［J］ International Archives of Photogrammetry，Remote Sensing and SpatialInformation Sciences，2007，36:413-418.

［10］ RONCAT A，WAGNER W，MELZER T，et al.Echo Detection and Localization in Fullwaveform Airborne Laser Scanner Data Using the Averaged Square Difference Function Estimator［J］.The Photogrammetric Journal of Finnland，2008，21（1）:62-75.

［11］ JUTZI B，STILLA U.Range Determination with Waveform Recording Laser Systems Using a Wiener Filter［J］.ISPRS JPRS，2006，61（2）:95-107.

［12］ WAGNER W，ULLRICH A，MELZER T，et al.From Single-Pulse to Full-Waveform Airborne Laser Scanners:Potential and Practical Challenges［C］∥Proceedings of International Society for Photogrammetry and Remote Sensing XXth Congress:XXXV.［S.l.］:ISPRS，2004.

［13］ DEMPSTER A P，LAIRD N M，RUBIN D B.Maximum Likelihood from Incomplete Data via the EM Algorithm［J］.J Royal Stat Soc Ser B-Methodol，1977，39:1-38.

［14］ 馬洪超，李奇.改進的EM模型及其在激光雷達全波形數據分解中的應用［J］.遙感學報，2009，13:35-41.

［15］ SUN G，RANSON K J.Modeling LiDAR Returns from Forest Canopies［J］.IEEE Trans Geosci Remote Sensing，2000，38:2617-2626.

［16］ BALTSAVIAS E P.A Comparison between Photogrammetry and Laser Scanning［J］.ISPRS Journal of Photogrammetry and Remote Sensing，1999，54（2/3）:83-94.