機(jī)載小光斑全波形LiDAR數(shù)據(jù)處理及應(yīng)用

周靜平，張愛武，王書民

（首都師范大學(xué)三維信息獲取與應(yīng)用教育部重點(diǎn)實(shí)驗(yàn)室，北京 100048）

一、引 言

機(jī)載激光雷達(dá)（airborne laser scanning，ALS）也指雷達(dá)或激光雷達(dá)，是一種通過激光器向地表發(fā)射紅外線短波脈沖，并由光電接收器接收背向散射的回波信號(hào)的主動(dòng)式光學(xué)遙感器［1］。作為一種先進(jìn)的新型對(duì)地觀測(cè)系統(tǒng)，它集激光、全球定位系統(tǒng)和慣性導(dǎo)航系統(tǒng)（inertial navigation system，INS）3項(xiàng)技術(shù)于一身，并用于快速獲取地面目標(biāo)的三維空間信息。

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

二、全波形LiDAR理論

1.激光傳感器

激光雷達(dá)的工作原理與傳統(tǒng)雷達(dá)遙感大體相同。激光測(cè)距系統(tǒng)是其核心組成部分，以激光波束掃描的方式測(cè)量激光器到目標(biāo)方向激光照射點(diǎn)的距離，即通過測(cè)量地面點(diǎn)的回波脈沖與激光的發(fā)射脈沖之間的時(shí)間延遲得到激光器到地面點(diǎn)之間的距離［4］。測(cè)距原理公式為

式中，R為激光器到目標(biāo)物的距離；c為光速；t為激光脈沖從激光器到目標(biāo)物傳輸?shù)耐禃r(shí)間［5］。機(jī)載激光雷達(dá)系統(tǒng)的R和t是可以相互轉(zhuǎn)換的，具體使用哪個(gè)可根據(jù)實(shí)際需求進(jìn)行選擇。激光雷達(dá)的測(cè)量方式是十分精確的，其測(cè)量精度常常可達(dá)到毫米級(jí)。再結(jié)合GPS記錄的坐標(biāo)信息和INS記錄的姿態(tài)角信息，通過聯(lián)合解算便可得到激光點(diǎn)的大地坐標(biāo)。

2.全波形激光雷達(dá)

全波形LiDAR即全數(shù)字波形激光雷達(dá)，采用數(shù)字化方式，不僅記錄若干次離散回波信號(hào)，而且將發(fā)射信號(hào)和回波信號(hào)均以很小的采樣間隔進(jìn)行采樣并記錄［6-8］。激光雷達(dá)記錄的回波波形是對(duì)光斑內(nèi)各個(gè)點(diǎn)的反射信號(hào)按時(shí)間先后順序的記錄，回波波形代表了激光后向散射的能量，也就是說，回波波形可以看做是回波強(qiáng)度信息在接收時(shí)間軸上的一個(gè)函數(shù)。

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

3.雷達(dá)方程

結(jié)合傳統(tǒng)雷達(dá)方程，激光雷達(dá)方程可表示如下［1］

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

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

式中，ti為往返的時(shí)間；為目標(biāo)i的振幅；Sp，i為回波脈沖的標(biāo)準(zhǔn)差。

式中，Ss是發(fā)射脈沖的標(biāo)準(zhǔn)差；為發(fā)射脈沖的振幅；Si是目標(biāo)i回波的標(biāo)準(zhǔn)差；σi為目標(biāo)i總的后向散射截面，是目標(biāo)地物不同散射截面的積分。目標(biāo)的散射截面與地表的粗糙度有關(guān)，可根據(jù)發(fā)射脈沖信號(hào)與不同目標(biāo)的后向散射截面的積分得到一個(gè)接收的波形。

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

三、波形分解方法和步驟

1.分解方法

波形分析技術(shù)最先由美國宇航局（NASA）用于其機(jī)載大光斑LiDAR系統(tǒng)上，激光雷達(dá)波形數(shù)據(jù)分析也是以美國為首進(jìn)行發(fā)展的。目前，對(duì)機(jī)載激光全波形數(shù)據(jù)進(jìn)行分解，國外研究得較早也較為深入，國內(nèi)研究最近幾年也逐漸開始。

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

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

本文采用基于高斯數(shù)學(xué)模型的波形分解法，對(duì)研究區(qū)波形數(shù)據(jù)進(jìn)行波形分解和擬合，采用L-M算法［1］并利用最小二乘法原理不斷迭代優(yōu)化，最終得到最優(yōu)的參數(shù)，即每個(gè)單個(gè)波形的振幅、距離、波寬。

2.分解步驟

1）數(shù)據(jù)預(yù)處理（波形數(shù)據(jù)讀取、噪聲去除、波形平滑）。

2）參數(shù)初始化（峰值點(diǎn)檢測(cè)與初始參數(shù)估計(jì)）。

3）參數(shù)提?。ㄒ罁?jù)擬合函數(shù)和分解算法）。

四、試驗(yàn)成果和分析

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

圖1 波形顯示界面

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

表1 回波數(shù)目分布情況

在Matlab平臺(tái)上，對(duì)研究區(qū)全波形數(shù)據(jù)進(jìn)行了顯示，并且利用激光點(diǎn)云數(shù)據(jù)處理專業(yè)軟件Terrasolid對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行了大致分類（見表2）。由于研究區(qū)域較大，為了方便顯示，只截取了其中一部分，效果如圖2所示。

表2 分類情況

從表1和圖1可以明顯看出，研究區(qū)的回波次數(shù)最多為7次，最少為1次。大部分為1次回波占87.814%；2次以上回波僅占12.186%。這是因?yàn)檠芯繀^(qū)位于城區(qū)，雖然地物類型較為豐富，但因建筑用地（房屋建筑、道路等）面積占據(jù)了大部分（占86.92%），這些地區(qū)較為平坦，激光對(duì)其無穿透；此研究區(qū)有也較多高大樹木、低矮灌木等（占13.08%），林木等地物具有葉間、枝間縫隙，激光可穿透。對(duì)比地物分類圖，可以明顯看出，建筑用地及道路等非植被區(qū)域和一次回波區(qū)大致重合，僅有1次回波；而有林木分布的植被區(qū)域和多次回波區(qū)大致吻合，回波為多次。整體上，也可以從表1和表2中1次回波的百分比和非植被區(qū)的百分比大致相同，多次回波的總百分比和植被區(qū)的百分比大致相等來進(jìn)一步得到說明。值得注意的是，研究區(qū)還有部分回波為0，這是因?yàn)榈匚锓瓷洳ǖ姆较虿辉诩す饨邮掌鞯囊曈騼?nèi)，導(dǎo)致接收不到地物的回波信息，或者是地物的反射率太低所致。

圖2 回波次數(shù)及地物分類圖

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

經(jīng)過波形分解，提取出了回波脈沖的距離、振幅和波寬等參數(shù)信息，結(jié)合點(diǎn)云分類信息，所成圖件如圖3所示。

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

圖3 回波參數(shù)及分類圖

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

圖4 振幅分布圖

五、結(jié)論與展望

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

機(jī)載小光斑全波形數(shù)據(jù)具有能夠快速獲取地表物體的精確三維空間信息和很多的地物特征信息的優(yōu)勢(shì)，未來還將繼續(xù)在地形測(cè)繪、農(nóng)林業(yè)、數(shù)字城市建設(shè)等諸多領(lǐng)域得到越來越廣泛和深入的應(yīng)用。

［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ī)載激光雷達(dá)數(shù)據(jù)的誤差分析及校正［J］.遙感學(xué)報(bào)，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］ 馬洪超，李奇.改進(jìn)的EM模型及其在激光雷達(dá)全波形數(shù)據(jù)分解中的應(yīng)用［J］.遙感學(xué)報(bào)，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.