全波形LiDAR數(shù)據(jù)分解的可變分量高斯混合模型及RJMCMC算法

趙泉華，李紅瑩，李 玉

遼寧工程技術(shù)大學(xué)測繪與地理科學(xué)學(xué)院，遼寧 阜新123000

1 引 言

與傳統(tǒng)激光雷達(dá)（light detection and ranging，LiDAR）系統(tǒng)記錄點云數(shù)據(jù)不同，全波形LiDAR系統(tǒng)能夠以波形的形式記錄一定高程范圍內(nèi)不同高程點上的后向散射能量［1－2］。記錄的波形信息可以提供一定范圍（光斑）內(nèi)的地表高程信息，以及不同高程點上反射能量。根據(jù)不同高程點反射能量的大小，能夠檢測出激光光束傳播過程中目標(biāo)垂直方向（高程方向）的分布情況［3－4］。根據(jù)獲取數(shù)據(jù)時激光束照射在地面光斑的大小，全波形LiDAR系統(tǒng)可分為大光斑系統(tǒng)和小光斑系統(tǒng)。大光斑LiDAR系統(tǒng)的光斑直徑一般為8～70 m，脈沖發(fā)射頻率為40～500 Hz，光斑之間間距在170 m左右，光斑分布稀疏；小光斑LiDAR系統(tǒng)的光斑直徑一般為0．3～1 m，脈沖發(fā)射頻率為40～500 k Hz，光斑間的間距小于1 m，光斑分布密集［5］。目前，常用的配有全波形數(shù)字化儀的LiDAR系統(tǒng)包括NASA機載的LVIS（laser vegetation imaging sensor）系 統(tǒng)［6－7］、SLICER（scanning LiDAR imager of canopies by echo recover）系 統(tǒng)［8－9］和 星 載 的 GLAS（geoscience laser altimeter system）系 統(tǒng)［10－11］等 大 光 斑LiDAR 系統(tǒng)；Optech ALTM3100、Leica ALS－70、Top Eye Mark Ⅱ、Riegl LMS－Q560 以 及Falcon Ⅲ等小光斑 LiDAR 系統(tǒng)［12－14］。

全波形LiDAR數(shù)據(jù)記錄的是地物在不同高程點上后向散射能量，因此全波形LiDAR數(shù)據(jù)不僅蘊含地物的高程信息，而且能夠反映地物的散射特性。波形分解是獲取這些信息的重要手段。波形分解是指用有限一致分布統(tǒng)計模型（Gaussian、Nakagami、Burr等）的組合來擬合原始波形數(shù)據(jù)。

目前，已經(jīng)有許多學(xué)者對全波形LiDAR數(shù)據(jù)的分解進行了研究。文獻［15］通過雷達(dá)方程證明高斯函數(shù)可以擬合全波形LiDAR數(shù)據(jù)，為分解波形數(shù)據(jù)提供了理論依據(jù)。文獻［16］提出了基于高斯模型的波形分解方法。該方法雖然能夠?qū)⒃疾ㄐ螖?shù)據(jù)分解為若干個高斯分量，但容易陷于局部最優(yōu)解，無法得到高精度的擬合結(jié)果。文獻［17］針對波形數(shù)據(jù)的高斯特性，對原始波形數(shù)據(jù)進行了分解，該方法雖然提高了擬合精度，但是容易遺漏部分重疊的返回波，分解結(jié)果不能充分反映地物沿高程方向的信息。

通過高斯分解，分析高斯分量數(shù)、均值和標(biāo)準(zhǔn)差等參數(shù)，能夠得到地物高程信息和后向散射特性。上述文獻中高斯分量數(shù)均為給定值，然而由于受地形起伏復(fù)雜、回波信號噪聲以及地物多次漫反射等諸多因素影響，很難人為確定波形分解數(shù)，并且容易遺漏部分重疊的返回波，降低波形擬合精度。因此，為了自動確定波形分解數(shù)，獲得精確的參數(shù)估計值，本文提出了一種基于可變分量高斯混合模型和可逆跳轉(zhuǎn)馬爾科夫鏈蒙特卡洛（reversible jump Markov chain Monte Carlo，RJMCMC）的全波形LiDAR數(shù)據(jù)分解方法。

2 算法描述

給定的全波形LiDAR數(shù)據(jù)集x＝｛xi；i＝1，2，…，n｝，y＝｛yi；i＝1，2，…，n｝，其中，xi為分布在有效高程L內(nèi)的采樣點，yi為在采樣點xi處的反射能量，n為采樣點個數(shù)。由此，全波形LiDAR數(shù)據(jù)可進一步表示為：z＝｛（xi，yi）；i＝1，2，…，n｝。一束激光的波形數(shù)據(jù)為由采樣點x作為橫坐標(biāo)，能量值y作為縱坐標(biāo)形成的點集，如圖1所示。

圖1 波形數(shù)據(jù)Fig．1 Waveform data

在波形分解過程中，首先假設(shè)理想波形為高斯混合模型，用理想波形來擬合實際波形數(shù)據(jù)；定義能量函數(shù)以刻畫二者的擬合程度；通過能量函數(shù)定義吉布斯分布，表征全波形LiDAR數(shù)據(jù)的或然率；根據(jù)貝葉斯定理構(gòu)建刻畫波形分解的后驗概率模型；設(shè)計RJMCMC算法模擬該后驗概率模型，從而確定波形分解個數(shù)及各高斯分量參數(shù)，完成波形分解。圖2所示為該算法總體流程，其中，實線框為完成的功能，虛線框為完成功能需要的模型和函數(shù)。

圖2 波形分解流程Fig．2 Waveform

2．1 理想波形模型

設(shè)理想波形為高斯混合模型，表達(dá)為m個高斯函數(shù)分量加權(quán)和

式中，＾yi（i＝1，2，…，n）為能量yi的理想值；wk為第k個高斯分量的權(quán)重，滿足和σk（k＝1，2，…，m）分別代表第k個高斯分量的均值和標(biāo)準(zhǔn)差。

2．2 或然率模型

為了控制理想模型與波形數(shù)據(jù)擬合程度定義能量函數(shù)。設(shè)參數(shù)集Φ＝｛μ，σ，w，m｝，其中，μ＝｛μ1，μ2，…，μm｝；σ＝｛σ1，σ2，…，σm｝；w＝｛w1，w2，…，wm｝。定義控制yi和＾yi之間差異程度的能量函數(shù)為U（z｜Φ），由試驗分析能量函數(shù)用絕對誤差和殘差二次型的平方根的效果相同，但采用絕對誤差可以減少計算量，因此，用絕對誤差代替殘差二次型的平方根，即

由能量函數(shù)U（z｜Φ）定義表征波形數(shù)據(jù)或然率的吉布斯分布h（z｜Φ）為

式中，Z為歸一化常數(shù)。

2．3 后驗概率模型

為了分解波形，需要構(gòu)建已知波形數(shù)據(jù)z條件下高斯分量的位置μ、標(biāo)準(zhǔn)差σ、權(quán)重w、個數(shù)m的聯(lián)合后驗概率密度函數(shù)。根據(jù)Bayes定理，已知波形z條件下的后驗概率為

式中的先驗概率討論如下：

（1）高斯分量位置μ。高斯分量的均值μk為對應(yīng)第k個分解波形最大反射能量的高程，稱為該高斯分量的位置。將m個高斯分量的位置按照增序排列為0＜μ1＜μ2＜…＜μm＜L。假設(shè)μk（k＝1，2，…，m）服從區(qū)間［0，L］上的均勻分布且相互獨立，則其聯(lián)合概率密度函數(shù)為

（2）高斯分量標(biāo)準(zhǔn)差σ。假設(shè)各高斯分量的標(biāo)準(zhǔn)差σk（k＝1，2，…，m）服從均值和標(biāo)準(zhǔn)差分別為μσ和σσ的正態(tài)分布且相互獨立，則其聯(lián)合概率密度函數(shù)為

（3）高斯分量權(quán)重w。假設(shè)各高斯分量的權(quán)重wk（k＝1，2，…，m）服從［0，1］上均勻分布并相互獨立，則其概率密度函數(shù)為

（4）高斯分量個數(shù)m。假設(shè)高斯分量的個數(shù)m服從均值為λ的泊松分布，則其概率密度函數(shù)為，其中，2≤m≤kmax，kmax為預(yù)先指定的最大高斯分量數(shù)。

綜上，式（4）的后驗概率可重寫為

2．4 波形分解

波形分解是通過對構(gòu)建的后驗概率模型進行模擬求解，本文采用RJMCMC算法進行后驗概率模擬。

2．4．1 RJMCMC算法

RJMCMC算法實現(xiàn)參數(shù)空間維數(shù)改變時，對參數(shù)集Φ的相關(guān)采樣滿足概率分布（式（5））。提取候選參數(shù)矢量Φ＊，根據(jù)文獻［18］，其接受率表示為

式中，u是為確保維度平衡而提出的隨機矢量，如狀態(tài)空間維數(shù)增加時，則｜Φ＊｜＝｜Φ｜＋｜u｜；rm（·）為 選 擇 移 動 類型 m 的 概 率；式 中 的Jacobian項，由參數(shù)集的改變，即由（Φ，u）→Φ＊引起的。

當(dāng)參數(shù)空間維數(shù)不變時，稱為 M－H（Metropolis－Hastings）算法，式（6）接受率簡化為

式中，Φl為總參數(shù)矢量Φ 的第l個元素，候選參數(shù)依其相應(yīng)的概率分布抽取。

2．4．2 基于RJMCMC方法的波形分解

采用RJMCMC算法對式（5）后驗概率模型進行模擬，為了實現(xiàn)完備采樣，設(shè)計4種移動操作類型：①更新高斯分量參數(shù)；②更新高斯分量權(quán)重w；③合并或分裂高斯分量；④生成或刪除高斯分量。

2．4．2．1 更新高斯分量參數(shù)

更新高斯分量參數(shù)包括更新高斯分量位置μ和更新高斯分量標(biāo)準(zhǔn)差σ，該操作中參數(shù)空間維數(shù)沒有發(fā)生改變，其過程包括：在m個高斯分量中隨機抽取一個高斯分量k∈｛1，2，…，m｝，該操作改變該高斯分量位置μk或標(biāo)準(zhǔn)差σk，而保持其他參數(shù)不變。候選和分別服從以當(dāng)前μk和σk為均值，σt和σs為標(biāo)準(zhǔn)差的正態(tài)分布。其中σt和σs分別為μk和σk的迭代步長，是預(yù)先定義的常數(shù)。當(dāng)σt過大時可以很快遍歷整個μk取值范圍，但是得到μk的精度較低；當(dāng)σt過小時能夠取到較高精度的μk，但是需要耗費較長時間。同理，σs具有相同的作用。通過試驗，當(dāng)σt和σs分別取50和1時得到的結(jié)果比較理想。則改變μ和σ 的接受概率分別為：a（μk）＝min｛1，Rμ｝，

2．4．2．2 更新高斯分量權(quán)重w

在更新高斯分量權(quán)重w的操作中，參數(shù)空間維數(shù)沒有發(fā)生改變，其過程包括：在m個高斯分量中隨機抽取兩個高斯分量，假設(shè)選取第k和k＋1個高斯分量，該操作改變wk和wk＋1，保持其他參數(shù)不變。隨機生成s∈（0，1），則改變高斯分量權(quán)重前后需滿足平衡條件＝（wk＋wk＋1）s和＝（wk＋wk＋1）（1－s）。則改變w的接受概率為：a（wk，）＝min｛1，Rw｝，其中

2．4．2．3 合并或分裂高斯分量

在合并或分裂高斯分量操作中，高斯分量的個數(shù)m發(fā)生變化，同時參數(shù)空間維數(shù)也發(fā)生改變。在分裂處理中，首先以等概率選擇高斯分量k，將其分裂成兩個高斯分量，分別記為k＊和（k＋1）＊，它們的位置、標(biāo)準(zhǔn)差以及權(quán)重分別為（μ＊k，將μk＋1、μk＋2、…、μm重新標(biāo)號為，σ、w作同樣的操作，則參數(shù)μ＝｛μ1，μ2，…，μm｝、σ＝｛σ1，σ2，…，σm｝和w＝｛w1，w2，…，wm｝分別變?yōu)棣蹋剑?，分裂后，高斯分量的位置均服從以當(dāng)前μk為均值，σt為標(biāo)準(zhǔn)差的正態(tài)分布，其中σt為預(yù)先定義的常數(shù)；標(biāo)準(zhǔn)差均服從以當(dāng)前σk為均值，σs為標(biāo)準(zhǔn)差的正態(tài)分布，其中σs為預(yù)先定義的常數(shù)；權(quán)重滿足平衡條件＝wk（1－u1），u1為保證維度平衡而設(shè)置的隨機數(shù)，u1∈（0，1）。則分裂高斯分量的接受概率為as＝min｛1，Rs｝，其中

式中，Jacobian項為1；rs和rm分別為分解高斯分量和合并高斯分量的總概率，rs＝sm／2m，rm＝mm／2（m＋1），其計算方法sm和mm分別為隨機選擇合并或分裂高斯分量的概率，當(dāng)m＝1時，m1＝0，s1＝1；當(dāng) m＝mmax時，smmax＝0、mmmax＝1（mmax是 m 的最大值）；當(dāng) m＝2，3，…，mmax－1時，sm＝mm＝0．5。分解高斯分量的過程包括：在m狀態(tài)下選擇分裂高斯分量，其概率為sm；然后在原有的m個高斯分量中以等概率隨機抽取一個高斯分量，如k，其概率為1／m；將對應(yīng)的第k個高斯分量以等概率（1／2）分離為第k＊個高斯分量或第（k＋1）＊個高斯分量。

合并高斯分量是分裂高斯分量的對偶操作。其過程包括：在m＋1狀態(tài)下選擇合并高斯分量操作，其概率為mm；然后在原有的m＋1個高斯分量中以等概率隨機抽取一個高斯分量，其概率為1／m＋1，然后在與其相鄰的兩個高斯分量中隨機抽取一個與其合并，其概率為1／2。由于合并高斯分量操作是分裂高斯分量操作的對偶操作，因此其接受概率為：ac＝min｛1，1／Rs｝。

2．4．2．4 生成或刪除高斯分量

在生成或刪除高斯分量操作中，高斯分量數(shù)m發(fā)生變化，同時參數(shù)空間維數(shù)也發(fā)生改變。對于生成高斯分量操作，首先以等概率選擇一個高斯分量k，設(shè)新生成的高斯分量的位置、標(biāo)準(zhǔn)差和權(quán)重分別為服 從 （μk，μk＋1）之間的均勻分布由其服從的先驗分布中抽取，權(quán)重滿足平衡條件＝wku2和為保證維度平衡而設(shè)置的隨機數(shù)，u2∈（0，1），則生成一個高斯分量的接受概率為ab＝min｛1，Rb｝，其中

式中，Jacobian項為1；rb＝bm／m、rd＝dm／m＋1，rb和rd分別為生成高斯分量和分解高斯分量的總概率，其計算方法為bm和dm的選取和前面的sm和cm相同。生成高斯分量的過程包括在m狀態(tài)下選擇生成操作，其概率為bm；在原有的m個高斯分量中以等概率隨機抽取一個高斯分量k，設(shè)生成的高斯分量標(biāo)號為（k＋1）＊，然后將μk＋1、μk＋2、…、μm重新標(biāo)號為其概率為1／m。

刪除高斯分量是生成高斯分量的對偶操作，其過程包括：在m＋1狀態(tài)下選擇刪除操作，其概率為dm；在原有m＋1個高斯分量中以等概率隨機抽取一個高斯分量k，并刪除該高斯分量，然后將位置μk＋1、μk＋2、…、μm＋1重新標(biāo)號為、…，其概率為1／m＋1。由于刪除高斯分量操作是生成高斯分量操作的對偶操作，因此其接受概率為ad＝min｛1，1／Rb｝。

3 試驗與分析

3．1 試驗數(shù)據(jù)

試驗數(shù)據(jù)采用ICESat－GLAS測高數(shù)據(jù)，來源于美國冰雪中心（national snow and ice data center，NSIDC）（http：∥nsidc．org／data／icesat／）獲取，數(shù)據(jù)采集時間為2003年10月05日至2005年10月23日。GLAS系統(tǒng)為首個極地軌道星載大光斑LiDAR系統(tǒng)，每秒發(fā)射40次激光脈沖，每個激光脈沖在地面上形成直徑約為70 m的光斑，光斑之間的間隔為170 m。每個光斑對應(yīng)一個波形，每個波形有544個采樣點。每個數(shù)據(jù)采集時間段對應(yīng)一個操作周期。L1A和L2A操作周期的波形數(shù)據(jù)采樣點間隔為1 ns（15 cm），有效高程范圍為81．6 m。從L2B操作周期以后，對波形進行了壓縮處理，前151個采樣點間隔增加到4 ns（60 cm），有效高程范圍增加到150 m［19－20］。ICESat－GLAS共提供了15種標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品（GLA01—GLA15），用戶可以根據(jù)不同需求和目的選擇相應(yīng)的標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品。試驗數(shù)據(jù)采用33版本的GLA01和GLA14兩種標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品。其中GLA01中記錄了光斑號、時間、脈沖強度、經(jīng)緯度等信息，GLA14中記錄高程、大地經(jīng)緯度、回波振幅等地表高度信息［21］。

本文的研究區(qū)針對城市區(qū)域，圖3所示為研究區(qū)的遙感影像和ICESat－GLAS光斑分布，圖3（a）中由上至下4條線分別為2004年5月20日、2005年5月22日、2005年10月23日和2003年10月05日的衛(wèi)星軌跡。由遙感圖像可以看出，該城區(qū)主要地物為道路、建筑物、廣場、公園、礦區(qū)等。本文選擇了小區(qū)、廣場、公園、郊區(qū)、礦區(qū)等幾個具有代表性的區(qū)域進行研究，圖3（b）—（f）分別為選擇5個典型區(qū)域的遙感影像，圓圈內(nèi)為光斑所覆蓋的范圍。圖4為5個研究區(qū)域?qū)?yīng)的ICESat－GLAS原始波形數(shù)據(jù)，表1為5個研究區(qū)域的數(shù)據(jù)信息。

圖3 研究區(qū)遙感影像及ICESat－GLAS光斑的分布Fig．3 Remote sensing image of study area and distribution for ICESat－GLAS shots

圖4 GLAS原始波形Fig．4 GLAS raw waveform

表1 數(shù)據(jù)信息Tab．1 Data information

3．2 數(shù)據(jù)預(yù)處理

為便于對不同返回波能量進行比較和加權(quán)，增強不同波形之間的可比性，需對波形數(shù)據(jù)進行歸一化處理，使能量值落入到［0，1］區(qū)間。

為了提高波形擬合的精度，需要對波形數(shù)據(jù)進行去噪處理。定義噪聲閾值為單個回波的4倍標(biāo)準(zhǔn)差，即：ε＝4σr，其中，σr為前100個高程點能量值的標(biāo)準(zhǔn)差。當(dāng)脈沖能量小于閾值時認(rèn)為是噪聲，將其能量視為0。

3．3 試驗結(jié)果及討論

首先，采用本文提出的算法對5個區(qū)域的波形數(shù)據(jù)進行波形分解試驗。圖5為最終分解結(jié)果。

圖5 分解結(jié)果Fig．5 Decomposition results

圖5（a）為區(qū)域1波形分解結(jié)果，分解成5個高斯分量。根據(jù)該區(qū)域遙感影像（圖3（b））可以看出該區(qū)域主要地物有高低不齊的平頂建筑物，車棚以及樹木。為了驗證地物覆蓋，對實地進行了勘測，并測得其平均高程分別為20、14．5、3和1 m。因此該地區(qū)分解的5個高斯分量分別代表高建筑物返回回波、低建筑物返回回波、車棚返回回波、樹木返回回波以及地面返回回波。試驗得到的地物高程分別為20．3、14．4、2．7和1．2 m。圖5（b）為區(qū)域2波形分解結(jié)果，分解成兩個高斯分量，根據(jù)該區(qū)域遙感影像（圖3（c））可以看出該區(qū)域主要地物為低矮植被，實測高程為0．3 m。因此該區(qū)域分解的兩個高斯分量分別為低矮植被返回波和地面返回回波。試驗得到的植被高程為0．3 m。圖5（c）為區(qū)域3波形分解結(jié)果，分解成7個高斯分量，根據(jù)該區(qū)域遙感影像（圖3（d））可以看出該區(qū)域垂直方向上主要地物有亭子和較高的植被，實測高程分別為8和4．5 m。因此該地區(qū)分解的7個高斯分量分別代表亭子返回回波，中間的幾個高斯分量為植被返回回波，最后一個高斯分量為地面返回回波。試驗得到亭子和植被的高程分別為8和4．4 m。圖5（d）為區(qū)域4波形分解結(jié)果，分解成兩個高斯分量，根據(jù)該區(qū)域遙感影像（圖3（e））可以看出該區(qū)域主要地物為尖頂?shù)牡徒ㄖ?，?jīng)實地勘測，其高程為3 m。因此該區(qū)域分解的兩個高斯分量分別為建筑物頂部返回回波和地面返回回波，試驗得到建筑物的高程為2．7 m。圖5（e）為區(qū)域5波形分解結(jié)果，分解成7個高斯分量，根據(jù)該區(qū)域遙感影像（圖3（f））可以看出該區(qū)域垂直方向上主要地物是地面傾斜的礦區(qū)，實測高差為8．5 m。因此該地區(qū)分解的7個高斯分量分別代表礦區(qū)不同高程地面返回回波。試驗得到地面傾斜礦區(qū)的高差為8．6 m。

通過分解后波形可以看出，低矮植被返回信號能量比較集中，振幅較高，脈沖較窄（圖5（b））；較高植被返回信號能量較分散，振幅較低，脈沖較寬（圖5（c））；平頂建筑物返回信號能量較大，振幅較高，脈沖較窄（圖5（a））；尖頂建筑返回信號能量較小，振幅較低，脈沖較寬（圖5（d）；地面傾斜的礦區(qū)根據(jù)地面傾斜程度的不同返回信號能量、振幅、脈沖的寬度都有所不同（圖5（e）。

從圖5的分解結(jié)果可以看出，對于不同地物類型區(qū)域的波形數(shù)據(jù)，由本文算法得到的理想波形都能很好地擬合原始波形數(shù)據(jù)，并且分解的波形能夠與實際地物相對應(yīng)，能夠得到地物沿高程方向的散射特性，說明本文算法具有可行性和有效性。

圖6 10 000次迭代誤差變化曲線Fig．6 Change in error during 10 000 iterations

試驗進行了10 000次迭代，在迭代循環(huán)過程中，誤差很快減小。圖6為5個研究區(qū)域10 000次迭代誤差變化曲線。由圖6可以發(fā)現(xiàn)，前4000次迭代誤差有一定波動，在后6000次迭代中誤差達(dá)到穩(wěn)態(tài)，為了使試驗結(jié)果更加準(zhǔn)確，忽略前4000次的迭代結(jié)果，比較后6000次高斯分量個數(shù)，并計算其分別出現(xiàn)的概率（表2為區(qū)域1不同高斯分量個數(shù)出現(xiàn)的概率），將出現(xiàn)概率最大的高斯分量個數(shù)作為最終分解的波形個數(shù)。

表2 區(qū)域1不同高斯分量個數(shù)出現(xiàn)的概率Tab．2 Probability of different number of Gaussian components for area 1

該算法基于 Matlab 7．12．0開發(fā)，運行環(huán)境為Intel Core i5－3470－3．20 GHz－4 G，Window7操作系統(tǒng)。以區(qū)域1為例，該區(qū)域有544個采樣點，一次迭代完成5種移動操作的時間約為0．02 s。因此，完成10 000次迭代總共的操作時間約為4 min。更新高斯分量位置、更新高斯分量標(biāo)準(zhǔn)差、更新高斯分量權(quán)重、合并或分裂高斯分量以及生成或刪除高斯分量分別占用時間的比例約為：10%、6%、4%、45%以及35%。其中，合并或分裂高斯分量和生成或刪除高斯分量兩種操作中不僅需要判斷是否滿足進行該操作的條件，而且還要判斷操作后高斯分量的均值、標(biāo)準(zhǔn)差以及權(quán)重是否滿足條件，因此，需要的操作時間略長。

為了驗證本文算法的可靠性和準(zhǔn)確性，將本文方法與Matlab提供的CFTOOL高斯擬合方法進行比較。由于CFTOOL高斯擬合方法不能自動確定分解成高斯分量數(shù)，所以將本文算法得到的高斯分量數(shù)作為該方法高斯分量數(shù)，然后將原始波形數(shù)據(jù)分解為固定數(shù)的高斯分量之和。圖7為CFTOOL高斯擬合方法的波形分解結(jié)果。表3分別列出了兩種方法得到的5個試驗區(qū)域高斯分布的均值估計值和標(biāo)準(zhǔn)差估計值。對比兩種分解方法結(jié)果，本文方法不僅能自動獲得高斯分量數(shù)，還可以較好地將波形分解，分解的波形能夠與不同高程地物相對應(yīng)，同時得到的高斯曲線能夠較好地擬合原始波形。CFTOOL高斯擬合方法不能自動獲得高斯分量數(shù)，雖然分解波形能較好地擬合原始波形，但未能分解出與不同高程地物對應(yīng)的波形分量，如區(qū)域1的第3個高斯分量標(biāo)準(zhǔn)差估計和區(qū)域2的第2個高斯分量標(biāo)準(zhǔn)差估計值過大，振幅接近于平值，并且對部分波形未能擬合，特別是有細(xì)微能量變化的波形，如區(qū)域3的第1個脈沖和區(qū)域4的第1個脈沖，區(qū)域5的第2個脈沖。因此，CFTOOL高斯擬合方法分解的波形不能與地物相對應(yīng)，無法反映地物沿高程方向的散射特性，其原因在于CFTOOL高斯擬合方法沒有對分解波形的參數(shù)進行約束，因此出現(xiàn)標(biāo)準(zhǔn)差估計值較大等情況。

圖7 CFTOOL方法分解結(jié)果Fig．7 Decomposition results from CFTOOL

表3 參數(shù)估計值Tab．3 Parameter estimation

為了對波形分解結(jié)果進行定量評價，分別計算每個試驗區(qū)域原始波形與擬合波形之間的相關(guān)系數(shù)ρ和KS距離［22］，計算公式分別為

式中，yi為原始波形中第i個高程點的能量值為擬合波形中第i個高程點的能量值和為其分別的均值。相關(guān)系數(shù)ρ∈［－1，1］，是衡量兩個隨機變量之間線性相關(guān)程度的指標(biāo)，相關(guān)系數(shù)越大，說明擬合數(shù)據(jù)與原始數(shù)據(jù)相關(guān)性越強，波形擬合的精度越高。KS∈［0，1］，KS距離能夠檢測出脈沖的丟失以及擬合數(shù)據(jù)與原始數(shù)據(jù)之間局部的變化。KS距離越小，說明波形擬合的精度越高。當(dāng)分解的波形能夠準(zhǔn)確地擬合原始波形數(shù)據(jù)時ρ＝1，KS＝0。

表4為計算結(jié)果，根據(jù)計算結(jié)果可以看出，得到的相關(guān)系數(shù)均在0．98以上，KS距離均在0．2以下，說明本文算法具有有效性和準(zhǔn)確性。

表4 相關(guān)系數(shù)及KS距離Tab．4 Cross－correlation coefficient and KS distance

4 結(jié) 論

本文提出了一種基于可變分量混合高斯模型和RJMCMC的全波形LiDAR數(shù)據(jù)分解方法，首先假定波形數(shù)據(jù)服從高斯分布，然后對LiDAR波形數(shù)據(jù)實現(xiàn)可變個數(shù)的分解。本文方法不僅能自動確定波形分解數(shù)，同時還可以有效地獲取最優(yōu)的模型參數(shù)估計。利用本文算法分別對不同地物類型的區(qū)域進行了波形分解試驗，并和Matlab中CFTOOL的高斯擬合方法進行了對比，定性和定量的測試結(jié)果驗證了本文方法的有效性、可靠性和準(zhǔn)確性。雖然本文算法能夠準(zhǔn)確地確定波形分解的個數(shù)，實現(xiàn)波形的分解，并且分解的結(jié)果能夠與地物相對應(yīng)，同時具有很高的擬合精度，但是對于一些非對稱、傾斜或者拖尾的波形，本文的方法還不能非常準(zhǔn)確地擬合。在未來的工作中，對如何進一步提高波形擬合精度，使算法更普遍地適用于各種形狀的波形，還需要繼續(xù)從理論上進行研究，如選擇其他的隨機模型來模擬原始波形數(shù)據(jù)，使分解的波形更準(zhǔn)確地擬合原始波形數(shù)據(jù)，從而提高波形擬合精度。

［1］ ZUO Zhiquan，ZHANG Zuxun，ZHANG Jianqing．A Highquality Filtering Method with Adaptive TIN Models for Urban LiDAR Points Based on Priori－knowledge［J］．Acta Geodaetica et Cartographica Sinica，2012，41（2）：246－251．（左志權(quán)，張祖勛，張劍清．知識引導(dǎo)下的城區(qū)Li－DAR點云高精度三角網(wǎng)漸進濾波方法［J］．測繪學(xué)報，2012，41（2）：246－251．）

［2］ SUI Lichun，ZHANG Yibin，LIU Yan，et al．Filtering of Airborne LiDAR Point Cloud Data Based on the Adaptive Mathematical Morphology ［J］． Acta Geodaetica et Cartographica Sinica，2010，39（4）：390－396．（隋立春，張熠斌，柳艷，等．基于改進的數(shù)學(xué)形態(tài)學(xué)算法的LiDAR點云數(shù)據(jù)濾波［J］．測繪學(xué)報，2010，39（4）：390－396．）

［3］ PARRISH C E，NOWAK R D．Improved Approach to LiDAR Airport Obstruction Surveying Using Full－waveform Data［J］．Journal of Surveying Engineering，2009，135（2）：72－78．

［4］ SUN G，RANSON K J．Modeling LiDAR Returns from Forest Canopies［J］．IEEE Transactions on Geoscience and Remote Sensing，2000，38（6）：2617－2626．

［5］ MALLET C，BRETAR F．Full－waveform Topographic LiDAR： State－of－the－art ［J］． ISPRS Journal of Photogrammetry and Remote Sensing，2009，64（1）：1－16．

［6］ BLAIR J B，RABINE D L，HOFTON M A．The Laser Vegetation Imaging Sensor：A Medium－altitude，Digitisationonly，Airborne Laser Altimeter for Mapping Vegetation and Topography［J］．ISPRS Journal of Photogrammetry and Remote Sensing，1999，54（2－3）：115－122．

［7］ SUN Guoqing，RANSON K J，ZHANG Zhongjun．Forest Vertical Parameters from LiDAR and Multi－angle Imaging Spectrometer Data［J］．Journal of Remote Sensing，2006，10（4）：523－530．（孫國清，RANSON K J，張鐘軍．利用激光雷達(dá)和多角度頻譜成像儀數(shù)據(jù)估測森林垂直參數(shù)［J］．遙感學(xué)報，2006，10（4）：523－530．）

［8］ PANG Yong，LI Zengyuan，CHEN Erxue，et al．LiDAR Remote Sensing Technology and Its Application in Forestry［J］．Scientia Silvae Sinicae，2005，41（3）：129－136．（龐勇，李增元，陳爾學(xué)，等．激光雷達(dá)技術(shù)及其在林業(yè)上的應(yīng)用［J］．林業(yè)科學(xué)，2005，41（3）：129－136．）

［9］ MEANS J E，ACKER S A，HARDING D J，et al．Use of Large－footprint Scanning Airborne LiDAR to Estimate Forest Stand Characteristics in the Western Cascades of Oregon［J］．Remote Sensing of Environment，1999，67（3）：298－308．

［10］ BRENNER A C，ZWALLY H J，BENTLEY C R，et al．ATBD V3．0 Derivation of Range and Range Distributions from Laser Pulse Waveform Analysis for Surface Elevations，Roughness，Slope，and Vegetation Heights［J／OL］．http：∥glas．wff．nasa．gov．2000．

［11］ SPAANS K H，RONSE A L A B．Reconstruction of an ICESat Full Waveform by Terrestrial and Airborne Laser Scanning［J］．Earth and Planetary Observation Letters，2007，3（20）：1－9．

［12］ LI Qi，MA Hongchao．The Study of Point－cloud Production Method Based on Waveform Laser Scanner Data［J］．Acta Geodaetica et Cartographica Sinica，2008，37（3）：349－354．（李奇，馬洪超．基于激光雷達(dá)波形數(shù)據(jù)的點云生產(chǎn)［J］．測繪學(xué)報，2008，37（3）：349－354．）

［13］ WANG Jinhu，LI Chuanrong，ZHOU Mei．Airborne Full－waveform LiDAR Data Processing and Application［J］．Foreign Electronic Measurement Technology，2012，31（6）：71－75，79．（王金虎，李傳榮，周梅．機載全波形激光雷達(dá)數(shù)據(jù)處理及其應(yīng)用［J］．國外電子測量技術(shù)，2012，31（6）：71－75，79．）

［14］ HUG C，ULLRICH A，GRIMM A．Litemapper－5600－A Waveform－digitizing LiDAR Terrain and Vegetation Mapping System［J］．International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2004，36（8／W2）：24－29．

［15］ WAGNER W，ULLRICH A，DUCIC V，et al．Gaussian Decomposition and Calibration of a Novel Small－footprint Full－waveform Digitising Airborne Laser Scanner［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2006，60（2）：100－112．

［16］ 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．

［17］ MA Hongchao，LI Qi．Modified EM Algorithm and Its Application to the Decomposition of Laser Scanning Waveform Data［J］．Journal of Remote Sensing，2009，13（1）：35－41．（馬洪超，李奇．改進的EM模型及其在激光雷達(dá)全波形數(shù)據(jù)分解中的應(yīng)用［J］．遙感學(xué)報，2009，13（1）：35－41．）

［18］ GREEN P J．Reversible Jump Markov Chain Monte Carlo Computation and Bayesian Model Determination［J］．Biometrika，1995，82（4）：711－732．

［19］ WEN Hanjiang，CHENG Pengfei．Introduction to Principle of ICESAT／GLAS Laser Altimetry and Its Applications［J］．Science of Surveying and Mapping，2005，30（5）：33－35．（文漢江，程鵬飛．ICESAT／GLAS激光測高原理及其應(yīng)用［J］．測繪科學(xué)，2005，30（5）：33－35．）

［20］ SUN G，RANSON K J，KIMES D S，et al．Forest Vertical Structure from GLAS：An Evaluation Using LVIS and SRTM Data［J］．Remote Sensing of Environment，2008，112（1）：107－117．

［21］ NELSON R，RANSON K J，SUN G，et al．Estimating Siberian Timber Volume Using MODIS and ICESat／GLAS［J］．Remote Sensing of Environment，2009，113（3）：691－701．

［22］ MALLET C，LAFARGE F，ROUX M，et al．A Marked Point Process for Modeling LiDAR Waveforms［J］．IEEE Transactions on Image Processing，2010，19（12）：3204－3221．