專利名稱：一種基于改進(jìn)的em算法的激光雷達(dá)波形數(shù)據(jù)分解的方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種對激光雷達(dá)波形數(shù)據(jù)處理的方法，特別是一種通過改進(jìn)EM 算法實(shí)現(xiàn)激光雷達(dá)波形數(shù)據(jù)分解的方法，屬于機(jī)載激光雷達(dá)技術(shù)領(lǐng)域。
背景技術(shù)：
機(jī)載激光掃描是一個迅速成長的地形測繪技術(shù)。機(jī)載激光掃描大多數(shù)都重復(fù) 對地球表面發(fā)射紅外脈沖，通過光學(xué)接收器檢測返回傳感器的一部分能量。由一 個計(jì)時器來測量脈沖在激光傳感器與地球表面之間的往返時間，通過這個往返時 間可以求得傳感器到地表的距離，加上經(jīng)過集成的P0S系統(tǒng)和大約2000萬 4000 萬像素的數(shù)碼相機(jī)，使得該技術(shù)在一些使用傳統(tǒng)攝影測量方法進(jìn)行地形測繪比較 困難的地區(qū)，如植被覆蓋區(qū)、海岸帶、島礁地區(qū)等等，發(fā)揮著不可替代的作用。
第一個商業(yè)機(jī)載激光掃描儀只能記錄一個后向散射脈沖的時間。如果在激光 光斑里只有一個目標(biāo)那僅僅一個脈沖記錄是足夠的，在這種情況下返回脈沖可以 簡單地說明地物。但是對小的激光光斑(0.2-2米)來說，在激光脈沖的傳播路 徑里可能有若干地物產(chǎn)生各自的后向散射脈沖。已經(jīng)有更加先進(jìn)的激光掃描儀能 夠記錄多次回波，但這種工作方式使用戶無法獲得任何與有關(guān)設(shè)備本身相關(guān)的一 些信息的，如如何根據(jù)回波信號定位，電子設(shè)備和地物結(jié)構(gòu)對獲取的回波信號的 形狀和大小有何影響，回波信號如何被量化成幾次離散脈沖信號等。另一方面， 機(jī)載激光雷達(dá)的光斑設(shè)計(jì)得越來越小，這樣每束激光產(chǎn)生多次回波的數(shù)量也進(jìn)一 步增加，這對脈沖探測和量化方法提出了更大的挑戰(zhàn)。曾有文獻(xiàn)指出，探測和量 化方法不同，可以引起成果的誤差。
解決上述問題的途徑是將發(fā)射信號和回波信號均以很小的采樣間隔進(jìn)行采 樣并記錄，而不僅僅是記錄若干次離散的回波信號。這樣的采樣記錄方式即所謂 的全波形數(shù)字化記錄，這種類型的激光雷達(dá)系統(tǒng)稱為waveform-digitizing LIDAR。事實(shí)上早在上世紀(jì)90年代，NASA開發(fā)的一些機(jī)載激光雷達(dá)系統(tǒng)己經(jīng)具 備這個能力，比如SLICER(Scanning Lidar Imager of Canopies by Echo Recover) 和LVIS (Laser Vegetation Imaging Sensor)。還有一些星載激光雷達(dá)系統(tǒng)如 GLAS也具有全波形記錄能力。但是與硬件發(fā)展相比，波形數(shù)據(jù)的分析和處理方法的研究卻相對滯后。這一方面是由于波形數(shù)據(jù)尚未真正普及，另一方面由于波
形數(shù)據(jù)的分析與具體應(yīng)用關(guān)系密切，分析方法要面向應(yīng)用。

發(fā)明內(nèi)容
本發(fā)明的目的就是提供一種處理激光雷達(dá)波形數(shù)據(jù)的方法，根據(jù)激光雷達(dá)波
形數(shù)據(jù)符合高斯分布的特征，利用改進(jìn)的EM (Expectation-Maximization)算
法高斯分解激光雷達(dá)波形數(shù)據(jù)，生產(chǎn)出高質(zhì)量的點(diǎn)云數(shù)據(jù)并得到森林參數(shù)。
實(shí)現(xiàn)本發(fā)明目的采用的技術(shù)方案 一種基于改進(jìn)EM算法的激光雷達(dá)波形數(shù)
據(jù)分解的方法，包括以下步驟
(1) 依據(jù)雷達(dá)等式的推導(dǎo)，用高斯函數(shù)模擬激光雷達(dá)的波形數(shù)據(jù)；
(2) 通過使用去噪閾值和平滑算法去除噪音；
(3) 在去噪后會出現(xiàn)毛刺現(xiàn)象，將毛刺作為噪音再進(jìn)行一次去噪處理；
(4) 在去除噪音后，通過波形維護(hù)來優(yōu)化波形；
(5) 對完全去除噪音后的激光雷達(dá)波形數(shù)據(jù)求梯度算子，把檢測到的若干激光 雷達(dá)波形極大值作為高斯函數(shù)的初始值；
(6) 通過EM算法對高斯函數(shù)參數(shù)做最大相似評估，將波形振幅當(dāng)作一個權(quán)值 加入到EM算法原始公式的分子和分母上；
(7) 用改進(jìn)的EM算法高斯分解波形數(shù)據(jù)，由于激光雷達(dá)波形數(shù)據(jù)會產(chǎn)生波形 重疊的情況，運(yùn)用最小距離法通過計(jì)算高斯函數(shù)期望值之間的差值來決定重疊波 形的高斯分解，所述改進(jìn)的EM算法即為將波形振幅當(dāng)作一個權(quán)值加入到EM算法 原始公式的分子和分母上；
(8) 提取分解激光雷達(dá)波形得到的高斯函數(shù)，高斯函數(shù)期望值為激光雷達(dá)波形 的位置，高斯函數(shù)的均方差為激光雷達(dá)波形的寬度；
(9) 通過應(yīng)用高斯函數(shù)期望值求解激光雷達(dá)點(diǎn)云數(shù)據(jù)三維坐標(biāo)，均方差求解森 林參數(shù)等地物屬性。
在步驟(2)和(3)中，去噪處理包括以下步驟-
(1) 通過求得發(fā)射脈沖波形的背景噪音平均值和均方差來確定返回脈沖波形的 去噪閾值，在返回脈沖波形里將小于去噪閾值的脈沖振幅賦值為零；
(2) 應(yīng)用平滑算法把波形曲線上的噪音抖動部分去掉。 在步驟(4)中，波形維護(hù)包括以下步驟
(1) 檢測去噪后每一段波形數(shù)據(jù)的位置；
(2) 在經(jīng)過平滑的原始波形數(shù)據(jù)上找到相應(yīng)的波形段，并把這些波形段兩端的 數(shù)據(jù)恢復(fù)到去噪后的波形數(shù)據(jù)里。
本發(fā)明通過用高斯分解的方式對激光雷達(dá)單個光斑范圍里的地物反射回波脈
沖進(jìn)行量化描述，實(shí)現(xiàn)對地物位置、反射率、粗糙度等一系列屬性提取，且利用
改進(jìn)的EM算法高斯分解激光雷達(dá)波形數(shù)據(jù)，生產(chǎn)出高質(zhì)量的點(diǎn)云數(shù)據(jù)并得到森
林參數(shù)。


下面結(jié)合附圖和實(shí)施例對本發(fā)明作進(jìn)一步說明。
圖l是原始波形數(shù)據(jù)。
圖2是帶有毛刺的去噪波形。
圖3是平滑后的去噪波形。
圖4是預(yù)處理后的波形。
圖5是維護(hù)前的波形。
圖6是維護(hù)后的波形。
圖7是波形的一階導(dǎo)數(shù)曲線。
圖8是設(shè)定初值后的波形圖。
圖9是基于EM算法波形高斯分解的初始流程圖。
圖10是基于EM算法的波形模擬結(jié)果。
圖11是改進(jìn)后的EM算法模擬的波形。
圖12是SLICER系統(tǒng)提取的地表位置。
圖13是改進(jìn)后的EM算法模擬的SLICER波形。
圖14是基于改進(jìn)后的EM算法波形高斯分解的最終流程圖。
具體實(shí)施例方式
本發(fā)明提供的基于改進(jìn)的EM算法的激光雷達(dá)波形數(shù)據(jù)處理的方法包括以下
步驟
(1)根據(jù)雷達(dá)等式的推導(dǎo)，用高斯函數(shù)模擬激光雷達(dá)的波形數(shù)據(jù)是一個最佳 方案，激光雷達(dá)原始波形數(shù)據(jù)如圖l所示，橫坐標(biāo)為波形數(shù)字化的采樣數(shù)，縱坐 標(biāo)為采樣的振幅值。
(2) 去除波形數(shù)據(jù)的噪音。波形在采集的過程中由于多方面的原因會存在噪 音，所以為一個抖動的曲線，圖1中振幅很小且抖動的部位即為噪音。先設(shè)定一 個閾值來去除噪音，去噪分兩步來處理第一步通過求得發(fā)射脈沖波形的背景噪 音平均值和均方差來確定返回脈沖波形的去噪閾值 。^ ，在返回脈沖波形里小
于去噪閾值 _的脈沖振幅賦值為零。由于波形的噪音振幅隨機(jī)性非常的大，
求取的去噪閾值a"。&無法把所有的噪音都去掉，造成去噪后的波形有毛刺，如
圖2所示。第二步應(yīng)用平滑算法把波形曲線上的噪音抖動部分去掉，如圖3所示， 去除噪音后的波形曲線是一條流暢平滑的信號曲線，這里可以看到有效波形部分 的噪音通過平滑算法被去掉，毛刺部分也因?yàn)槠交惴ǘ粔合氯ィ桃廊?沒有去除。
(3) 由于有些波形數(shù)據(jù)的噪音部分抖動得很厲害，在去噪后會出現(xiàn)很多毛刺 現(xiàn)象，檢測出毛刺，并將它們作為噪音再進(jìn)行一次去噪處理。再用去噪閾值0""^
做一次去噪處理，最終得到的結(jié)果如圖4所示。
(4) 如圖5所示，由于有些地物的反射率很低或者由于遮擋打在地物上的能 量很弱使其反射很小的能量，在去除噪音的時候使有效波形損失過大不利于后續(xù) 處理。對去噪后的波形進(jìn)行維護(hù)分兩步進(jìn)行第一步檢測去噪后每一段波形數(shù)據(jù) 的位置，第二步在經(jīng)過平滑的原始波形數(shù)據(jù)上找到這些相應(yīng)的波形段，并把這些 波形段兩端一定范圍的數(shù)據(jù)恢復(fù)到去噪后的波形數(shù)據(jù)里，圖6為維護(hù)后的波形。 圖5和圖6中的橫線為去噪閾值 。^ 。
(5) 使用EM算法之前，A.的初始值為波形的局部極大值，如圖6所示的位
置，因此可運(yùn)用梯度算子對波形求一階導(dǎo)數(shù)，圖7是圖6的一階導(dǎo)數(shù)曲線，左 邊大于0右邊小于0的采樣點(diǎn)就是局部極大值所在的位置，并由極大值的個數(shù)確
定初始總的高斯函數(shù)數(shù)量。A.的初始值是讓所有的J個高斯分支有相同的權(quán)， 在本實(shí)施例中C7y的初始值設(shè)置為1。如圖8所示，線1是預(yù)處理后的波形數(shù)據(jù)，
線2是設(shè)定初始值后的高斯函數(shù)曲線。
(6) 由于回波波形可以看作是若干高斯函數(shù)的疊加，波形的分解即為混合高 斯分布，通過EM算法完成高斯混合密度的參數(shù)估計(jì)。
EM算法計(jì)算巧、//;、 0";最優(yōu)值的原始公式為
<formula>complex formula see original document page 7</formula> (1)<formula>complex formula see original document page 7</formula> (2) <formula>complex formula see original document page 7</formula> (3)
<formula>complex formula see original document page 7</formula> (4)
其中Qu是采樣i屬于高斯函數(shù)分支j的權(quán)重，/;.(x)是高斯概率密度函數(shù)，
a.是/;.0c)的權(quán)，a.是高斯函數(shù)的期望值，。是高斯函數(shù)的標(biāo)準(zhǔn)差，通過求得
的//;、。.初始值得到初始高斯函數(shù)力(x) iv(^., 2)，波形由公式
/00 = 1>,.></;.(力來描述。通過迭代公式(i) (4)不斷地調(diào)節(jié)巧、^.、 .
這三個參數(shù)得到最優(yōu)值。波形數(shù)據(jù)分解的初始流程圖如圖9所示。得到波形分解 的結(jié)果圖如圖IO所示，左邊的高斯模擬波形是合理的，但右邊存在很大的偏差。
波形由一系列簡單的高斯分布組成，則混合分布的數(shù)學(xué)表示式(即采樣后得
到的波形)可以由下式表示
其中a是高斯函數(shù)的數(shù)量，力oc)是高斯概率密度函數(shù)，^是y;(力的權(quán)，表
示該分布在混合分布中占的比重，滿足0<巧<1，^巧=1， ^是高斯函數(shù)期
望值，^.是高斯函數(shù)的標(biāo)準(zhǔn)差。對每個部分J'，計(jì)算得到的a.表示回波的位置， 0";表示回波的寬度。參數(shù)a、 //;、 均可以由EM算法估計(jì)得到。先對數(shù)據(jù)進(jìn) 行預(yù)處理，并對公式(1) (4)的參數(shù)定義初始值。由于初始值是預(yù)測的，f (x) 也是依據(jù)初始值^、 得到的高斯函數(shù)，然后用f(x)來計(jì)算2^，但在后來的公式(1) - (4)中還是僅僅用預(yù)測的^來計(jì)算，而不考慮實(shí)際的振幅值來擬合
波形數(shù)據(jù)得不到一個理想結(jié)果。所以從物理意義上來說在公式(1) (4)里加 上一個振幅值m。以下為公式推導(dǎo)
將公式(2)代入到公式(3) (4)中得到<formula>complex formula see original document page 8</formula>
(5)
<formula>complex formula see original document page 8</formula>(6)
現(xiàn)在同時在分母和分子加入振幅A^， iV,相當(dāng)于一個權(quán)值，用來約束&、 a ^值，得到公式如下
<formula>complex formula see original document page 8</formula>(7)
<formula>complex formula see original document page 8</formula>(8)
把公式(7) - (8)還原為EM算法的形式
<formula>complex formula see original document page 8</formula>
n是波形中采樣的數(shù)量， 是第i次采樣的振幅。圖11為用改進(jìn)的EM算法得
到的結(jié)果，可以看出其很好地模擬了波形數(shù)據(jù)。
將改進(jìn)的EM算法得到的波形參數(shù)與SLICER系統(tǒng)得到的波形參數(shù)進(jìn)行比較SLICER提供了一個峰值提取算法，如圖12所示，三條豎線為系統(tǒng)提取的 Grstart、 Grpeak、 GrEnd三個峰值，分別代表地表返回的開始、峰值和末端。 在Grpeak、 GrEnd之間有一個波峰，這個波峰是真正的地表返回，這束激光產(chǎn)生 了兩次明顯的回波首先其打在了一個濃密的低矮植被上產(chǎn)生一個振幅較大的回 波，然后打在地表。由于大部分激光能量被植被反射，所以地表反射的回波較小。 又由于植被與地表距離較近，使兩次回波發(fā)生了部分重疊。由此可見系統(tǒng)提取的 地表返回位置不夠精確。
用改進(jìn)的EM算法得到的模擬波形如圖13所示，虛線為EM算法模擬的兩次
回波，兩條縱向直線標(biāo)明了兩次回波波峰的位置，//,=57、 a,=7.3、 //2=83， cr2=7.3。將圖13與圖12進(jìn)行對比，圖13得到的結(jié)果精度更高，正確性更好。 所以利用改進(jìn)的EM算法波形分解可以獲得更高的精度和更好的正確性。
(7) 在找到最優(yōu)值后，采用1^11而攀)=1^11(//7.+1-//J方法，確定得到的A:值
是否是最佳的，即用yt個高斯函數(shù)描述波形是否合適。圖14為是基于改進(jìn)后的 EM算法波形高斯分解的最終流程圖。
(8) 得到波形高斯分解的期望值后，依據(jù)波形數(shù)據(jù)起始點(diǎn)的坐標(biāo)可以計(jì)算點(diǎn) 云數(shù)據(jù)，波形分解得到的點(diǎn)云數(shù)據(jù)精度要比系統(tǒng)RIEGLLiteM即per5600的要高， 波形分解的點(diǎn)云數(shù)據(jù)也更加有層次感，奇異點(diǎn)的數(shù)量也明顯減少，其可以為DTM 生產(chǎn)、城市建模等提供優(yōu)質(zhì)的數(shù)據(jù)源。高斯分解得到的點(diǎn)云無論是數(shù)量還是質(zhì)量 都優(yōu)于系統(tǒng)RIEGLLiteMapper5600提取的點(diǎn)云數(shù)據(jù)。
(9) 在森林覆蓋地區(qū)，用高斯分解的最后一次回波到反射波形開始處之間的 距離來表示樹高；用激光打在樹冠上的脈沖總量來描述冠層蓄積量；單次回波用
高斯分解檢測，將檢測出的波形數(shù)據(jù)用公式地表反射率<formula>complex formula see original document page 9</formula>計(jì)算地
表反射率；將檢觀j出多次回波的數(shù)據(jù)用公式
植被反射率=<formula>complex formula see original document page 9</formula>計(jì)算植被反射率，
由于影響反射率的因數(shù)很多，用地表平均反射率來計(jì)算植被平均反射率會引起誤
差，但這個平均反射率在局部地區(qū)還是相當(dāng)有代表性的，其結(jié)果可以用于分類；
應(yīng)用公式森林簡度<formula>complex formula see original document page 9</formula>平均反射率來計(jì)算森林郁閉度。
權(quán)利要求
1.一種基于改進(jìn)的EM算法的激光雷達(dá)波形數(shù)據(jù)分解的方法，其特征在于包括以下步驟(1)依據(jù)雷達(dá)等式的推導(dǎo)，用高斯函數(shù)模擬激光雷達(dá)的波形數(shù)據(jù)；(2)通過使用去噪閾值和平滑算法去除噪音；(3)在去噪后會出現(xiàn)毛刺現(xiàn)象，將毛刺作為噪音再進(jìn)行一次去噪處理；(4)在去除噪音后，通過波形維護(hù)來優(yōu)化波形；(5)對完全去除噪音后的激光雷達(dá)波形數(shù)據(jù)求梯度算子，把檢測到的若干激光雷達(dá)波形極大值作為高斯函數(shù)的初始值；(6)通過EM算法對高斯函數(shù)參數(shù)做最大相似評估，將波形振幅當(dāng)作一個權(quán)值加入到EM算法原始公式的分子和分母上；(7)用改進(jìn)的EM算法高斯分解波形數(shù)據(jù)，由于激光雷達(dá)波形數(shù)據(jù)會產(chǎn)生波形重疊的情況，運(yùn)用最小距離法通過計(jì)算高斯函數(shù)期望值之間的差值來決定重疊波形的高斯分解，所述改進(jìn)的EM算法即為將波形振幅當(dāng)作一個權(quán)值加入到EM算法原始公式的分子和分母上；(8)提取分解激光雷達(dá)波形得到的高斯函數(shù)，高斯函數(shù)期望值為激光雷達(dá)波形的位置，高斯函數(shù)的均方差為激光雷達(dá)波形的寬度；(9)通過應(yīng)用高斯函數(shù)期望值求解激光雷達(dá)點(diǎn)云數(shù)據(jù)三維坐標(biāo)，均方差求解森林參數(shù)等地物屬性。
2. 根據(jù)權(quán)利要求1所述基于改進(jìn)的EM算法的激光雷達(dá)波形數(shù)據(jù)分解的方法，其 特征在于步驟(2)和(3)中，去噪處理包括以下步驟(2-l)通過求得發(fā)射脈沖波形的背景噪音平均值和均方差來確定返回脈沖波形 的去噪閾值，在返回脈沖波形里小于去噪閾值的脈沖振幅賦值為零； (2-2)應(yīng)用平滑算法將波形曲線上的噪音抖動部分去掉。
3. 根據(jù)權(quán)利要求1所述基于改進(jìn)的EM算法的激光雷達(dá)波形數(shù)據(jù)分解的方法，其特征在于步驟(4)中，波形維護(hù)包括以下步驟 (3-1)檢測去噪后每一段波形數(shù)據(jù)的位置；(3-2)在經(jīng)過平滑的原始波形數(shù)據(jù)上找到相應(yīng)的波形段，并將這些波形段兩端 的數(shù)據(jù)恢復(fù)到去噪后的波形數(shù)據(jù)里。
全文摘要
本發(fā)明公開了一種基于改進(jìn)EM算法的激光雷達(dá)波形數(shù)據(jù)分解的方法，屬于機(jī)載激光雷達(dá)技術(shù)領(lǐng)域。通過改進(jìn)的EM算法即將波形振幅當(dāng)作一個權(quán)值加入到EM算法原始公式的分子和分母上對激光雷達(dá)波形進(jìn)行分析，可以更加詳細(xì)地了解物體的縱向結(jié)構(gòu)，如表面傾斜、粗糙度、反射率。本發(fā)明在后處理過程中讓使用者自己提取三維坐標(biāo)獲得高精度的點(diǎn)云結(jié)果。本發(fā)明使用改進(jìn)的EM算法得到回波脈沖的位置和寬度，是一種性能可靠、精度較高的波形分解算法。
文檔編號G01S7/48GK101196562SQ20071016890
公開日2008年6月11日 申請日期2007年12月14日 優(yōu)先權(quán)日2007年12月14日
發(fā)明者奇 李, 馬洪超 申請人:武漢大學(xué)