鄭輯濤，張 濤

1．清華大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)系，北京100084；2．北京航空氣象研究所，北京100085

1 引 言

地形表面通常分布著建筑物、植被、車輛和人員等物體，激光掃描獲取地形點(diǎn)云后的基礎(chǔ)工作之一是把分布在地面上的點(diǎn)（地面點(diǎn)）與分步在其他目標(biāo)上的點(diǎn)（地物點(diǎn)）進(jìn)行區(qū)分。LiDAR點(diǎn)云濾波的目的就是區(qū)分地面點(diǎn)和地物點(diǎn)并濾除地物點(diǎn)得到地面點(diǎn)，這是建立高精度數(shù)字地面模型的關(guān)鍵環(huán)節(jié)。

近年來隨著三維激光掃描技術(shù)的進(jìn)步，很多學(xué)者對(duì)LiDAR點(diǎn)云濾波方法進(jìn)行了研究并提出了相關(guān)算法。根據(jù)文獻(xiàn)［1］中的總結(jié)，可以把現(xiàn)有的濾波方法大致分為以下幾類：

（1）基于坡度變化的方法。經(jīng)典的坡度變化濾波法［2］設(shè)計(jì)了兩點(diǎn)距離為自變量而最大高度差閾值為因變量的函數(shù)，根據(jù)兩點(diǎn)間實(shí)際高差與高差閾值比較區(qū)分地面點(diǎn)和地物點(diǎn)。

（2）三角網(wǎng)漸進(jìn)加密濾波法。文獻(xiàn)［3］首次提出了三角網(wǎng)漸進(jìn)加密濾波方法，其區(qū)分地物點(diǎn)和地面點(diǎn)的依據(jù)是當(dāng)前點(diǎn)與其垂足所在三角形之間的垂直距離和當(dāng)前點(diǎn)到該三角形3個(gè)頂點(diǎn)的連線與該三角形之間的夾角。三角網(wǎng)漸進(jìn)加密濾波法對(duì)森林地區(qū)的植被和城區(qū)的高大建筑物有較好的濾波效果，因此該方法得到較多關(guān)注［4－7］，并應(yīng)用于商業(yè)軟件TerraScan中。文獻(xiàn)［4］在三角網(wǎng)漸進(jìn)加密的基礎(chǔ)上引入了基于多核處理器的并行計(jì)算技術(shù)，提高了濾波效率以滿足快速成圖的應(yīng)用需求。該方法的缺點(diǎn)是濾波過程中初始網(wǎng)格選擇不當(dāng)容易積累誤差并且效率偏低，對(duì)于低矮植被和物體的濾除效果比較差。

（3）基于曲面擬合的方法。最有代表性的是文獻(xiàn)［8］中提出的最小二乘逼近法，該方法用最小二乘法擬合曲面，根據(jù)點(diǎn)到該曲面的距離分配權(quán)值，并根據(jù)閾值設(shè)計(jì)了一個(gè)權(quán)函數(shù)，能夠保證地面點(diǎn)權(quán)值為1而地物點(diǎn)權(quán)值為0，經(jīng)過多次迭代可以逐漸達(dá)到濾波的目的。此外也有研究者用Snake樣條［9］和 多 項(xiàng) 式［10－11］擬合 逼 近 地 形 表 面，此類方法在地形變化平緩的地區(qū)效果比較好，但對(duì)于大型建筑物往往會(huì)誤把建筑物頂面錯(cuò)判為地形表面。

（4）基于數(shù)學(xué)形態(tài)學(xué)的方法。文獻(xiàn)［12］最早提出了該方法，此后很多學(xué)者都在此基礎(chǔ)上提出了改進(jìn)的數(shù)學(xué)形態(tài)學(xué)濾波法［13－19］。該方法把圖像處理領(lǐng)域的數(shù)學(xué)形態(tài)學(xué)概念進(jìn)行擴(kuò)展，形成針對(duì)地形點(diǎn)云的腐蝕膨脹過程，構(gòu)造開運(yùn)算。其弊端是窗口的大小會(huì)直接影響濾波的效果，窗口過小導(dǎo)致濾波對(duì)大型建筑物失效，窗口過大會(huì)使地面上較大的凸起被錯(cuò)誤地濾除。文獻(xiàn)［13］在數(shù)學(xué)形態(tài)學(xué)開算子的基礎(chǔ)上，提出增加一個(gè)“帶寬”參數(shù)以提高其自適應(yīng)程度。此外對(duì)地形進(jìn)行腐蝕膨脹操作時(shí)整個(gè)窗口的高程取同一個(gè)值，使得地面細(xì)節(jié)丟失，從而造成地形失真。文獻(xiàn)［14］采取分步策略，在數(shù)學(xué)形態(tài)學(xué)粗濾波的基礎(chǔ)上，引入了Car模型對(duì)地形進(jìn)行了模型化處理，提高了濾波的穩(wěn)健性。

（5）其他方法。基于體素的濾波法［20］、利用多次回波特性的濾波方法［21］以及二次濾波法［22］等。

盡管濾波的方法已有多種，但試驗(yàn)對(duì)比表明［1，23］，準(zhǔn)確高效并且普遍適用的方法仍有待深入研究。

濾波過程會(huì)導(dǎo)致兩類錯(cuò)誤，一類是把本屬于地面的點(diǎn)錯(cuò)判為地物點(diǎn)濾除；另一類是把本屬于地物上的點(diǎn)錯(cuò)判為地面點(diǎn)保留。文獻(xiàn)［1］把這兩類錯(cuò)誤分別稱為Ⅰ型錯(cuò)誤和Ⅱ型錯(cuò)誤。通常地形表面都是連續(xù)緩和變化，少量的Ⅰ型錯(cuò)誤并不會(huì)導(dǎo)致地形的較大失真；而Ⅱ型錯(cuò)誤過多會(huì)引起地形的較大誤差。因此應(yīng)該在維持較小Ⅱ型錯(cuò)誤的基礎(chǔ)上盡量減少Ⅰ型錯(cuò)誤的發(fā)生。為了改善濾波精度、簡化濾波方法并提高濾波效率，本文給出一種改進(jìn)的濾波算法，其思路是：①沿同一方向等間距逐行掃描點(diǎn)云，獲取點(diǎn)序列構(gòu)成的掃描線，針對(duì)每條掃描線，采取半徑可變的圓環(huán)從面向地心一側(cè)滾過，滾動(dòng)過程中遠(yuǎn)離圓環(huán)邊緣的點(diǎn)通常是地物點(diǎn)可以濾除；②為便于采取B樣條擬合地形表面，對(duì)每條濾波后的掃描線等間距均勻采樣，得到規(guī)則分布的地表型值點(diǎn)，在此基礎(chǔ)上用均勻B樣條曲面擬合地形表面；③遍歷原始地形點(diǎn)云中的每一個(gè)點(diǎn)，計(jì)算其在擬合曲面上的投影高程，并與實(shí)際高程比較，根據(jù)兩者的差值區(qū)分地表點(diǎn)和地物點(diǎn)。技術(shù)流程如圖1所示。

圖1 算法流程Fig．1 Algorithm flow chart

2 掃描線濾波

2．1 掃描線獲取與掃描線濾波

通常情況下地面的變化比較平緩，獲取掃描線可以通過設(shè)定一組與地面垂直的等間隔平行平面，如圖2所示，然后給定閾值Δd，每一個(gè)平面附近到該平面的距離小于給定閾值Δd的點(diǎn)構(gòu)成該平面位置處的掃描線。閾值Δd可以取作點(diǎn)云平均間距的1／2。點(diǎn)云平均間距d可以通過式（1）計(jì)算得到

式中，n為點(diǎn)云中點(diǎn)的數(shù)量；L和W 分別表示點(diǎn)云所表示地形范圍的長和寬。

圖2 獲取掃描線Fig．2 Obtain scanning line

圖3所示為圖2所示地形點(diǎn)云中的一條掃描線，從該掃描線上可以清楚看出地面和植被、建筑等地物。

圖3 掃描線Fig．3 Scanning line

由于空洞的存在會(huì)導(dǎo)致掃描線上出現(xiàn)數(shù)據(jù)缺失，此時(shí)需要連接掃描線上空洞的兩端點(diǎn)，以點(diǎn)云平均間距為間隔進(jìn)行等間距插值，補(bǔ)充缺失的數(shù)據(jù)點(diǎn)。這使得掃描線濾波得以順利進(jìn)行。

掃描線濾波就是要濾除掃描線上的地物點(diǎn)，為此可以從面向地心的一側(cè)，用一個(gè)圓環(huán)，從掃描線的一端逐點(diǎn)依次滾向另一端，由于建筑物和植被等非地面點(diǎn)從地心一側(cè)觀察表現(xiàn)為向上凹陷，因此當(dāng)圓環(huán)半徑取值適當(dāng)時(shí)，圓環(huán)邊緣只能滾過地面點(diǎn)，未被圓環(huán)滾過的點(diǎn)被認(rèn)為是地物點(diǎn)并濾除。

2．2 圓心計(jì)算

對(duì)于任一條掃描線可以任選其一端為起始點(diǎn)，則另一端即為終點(diǎn)，如圖4所示，為了說明問題方便，不妨設(shè)掃描線的左端點(diǎn)為起始點(diǎn)，則當(dāng)圓環(huán)從面向地心的一側(cè)從左端起始點(diǎn)滾向右端終點(diǎn)時(shí)，圓環(huán)實(shí)際在作逆時(shí)針滾動(dòng)。

定義1 當(dāng)前點(diǎn)：由于掃描線上的點(diǎn)序列從起始點(diǎn)到終點(diǎn)有序排列，則圓環(huán)滾動(dòng)過程中任一時(shí)刻與圓環(huán)邊緣相接觸的點(diǎn)中最接近終點(diǎn)的一個(gè)稱為當(dāng)前點(diǎn)。

定義2 相鄰目標(biāo)點(diǎn)：掃描線上的點(diǎn)序列中，在圓環(huán)滾動(dòng)前進(jìn)方向上與當(dāng)前點(diǎn)相鄰的下一點(diǎn)稱為相鄰目標(biāo)點(diǎn)。

把起始點(diǎn)作為地表點(diǎn)和當(dāng)前點(diǎn)，圓環(huán)滾動(dòng)過程中依次判斷下一個(gè)相鄰目標(biāo)點(diǎn)是否被圓環(huán)的邊緣滾過，直至終點(diǎn)。滾過與否可以通過比較圓環(huán)半徑和點(diǎn)到圓心的距離來判斷。給定圓環(huán)半徑r，圓心的計(jì)算如圖4所示。

圖4 圓心計(jì)算Fig．4 Calculation circle center

假設(shè)當(dāng)前點(diǎn)為P1（x1，y1），滾向相鄰目標(biāo)點(diǎn)P2（x2，y2），可以計(jì)算矢量P1P2＝｛x2－x1，y2－y1｝，把該矢量順時(shí)針旋轉(zhuǎn)90°可得到過P1P2中點(diǎn)P0并指向圓心O的矢量n

因此得過P1P2中點(diǎn)P0并指向圓心O的參數(shù)方程

由于P0是P1P2的中點(diǎn)，因此有

根據(jù)

代入P0和圓心O得

代入P1P2的坐標(biāo)得

圓心O的對(duì)應(yīng)參數(shù)

代入?yún)?shù)方程可得圓心坐標(biāo)

當(dāng)從P1點(diǎn)滾向P2點(diǎn)時(shí)，首先根據(jù)式（9）計(jì)算圓心，然后判斷掃描線上是否存在到該圓心距離小于r的點(diǎn)：如果存在這樣的點(diǎn)，說明P2不能被圓環(huán)邊緣滾過，因此可以確定P2點(diǎn)是地物點(diǎn)，繼續(xù)取圓環(huán)前進(jìn)方向上與P2點(diǎn)相鄰的下一點(diǎn)作為相鄰目標(biāo)點(diǎn)，判斷能否從當(dāng)前點(diǎn)P1滾過；如果不存在這樣的點(diǎn)，說明圓環(huán)可以滾過P2點(diǎn)，則可以確定P2點(diǎn)是地面點(diǎn)，然后用P2更新P1作為當(dāng)前點(diǎn)，繼續(xù)判斷，直到圓環(huán)滾至終點(diǎn)。

2．3 可變半徑

對(duì)于實(shí)際地形，建筑物的尺寸大小不同，圓環(huán)半徑過大會(huì)把凸起的地形濾除，圓環(huán)半徑過小，會(huì)滾過大型建筑物，把部分建筑物表面誤判為地面。此外掃描線上相鄰點(diǎn)的間距并不均勻，如果采用固定半徑，當(dāng)相鄰點(diǎn)距大于圓環(huán)直徑時(shí)，圓環(huán)會(huì)漏過，表現(xiàn)在計(jì)算中是在計(jì)算圓心參數(shù)時(shí)出現(xiàn)復(fù)數(shù)開方的情況，導(dǎo)致計(jì)算中斷。因此為避免此類情況發(fā)生，可以設(shè)計(jì)可變半徑的圓環(huán)，半徑根據(jù)當(dāng)前點(diǎn)和相鄰目標(biāo)點(diǎn)之間的距離自動(dòng)調(diào)整?？勺儼霃降膱A環(huán)滾動(dòng)如圖5所示。

圖5 可變半徑Fig．5 Variable radius

為保持更多的地面細(xì)節(jié)，圓環(huán)的邊緣可以使之具有彈性變形，如同真實(shí)地面上具有彈性的輪胎能夠正常滾過微小的凹凸地形。

半徑的變化范圍可根據(jù)點(diǎn)云的平均間距和建筑物最大尺寸設(shè)定。不妨用r0表示點(diǎn)云的平均間距，作為圓環(huán)半徑的最小取值，用r1表示建筑物最大尺寸，作為圓環(huán)半徑的最大取值。當(dāng)相鄰的當(dāng)前點(diǎn)與相鄰目標(biāo)點(diǎn)之間的距離增大時(shí)，為避免出現(xiàn)計(jì)算中斷情況，半徑的變化應(yīng)該能夠較快的增大，然后逐漸穩(wěn)定在最大半徑附近，因此可以設(shè)計(jì)半徑計(jì)算公式（10），即

式中，d表示當(dāng)前點(diǎn)P1與相鄰目標(biāo)點(diǎn)P2之間的距離；Δ表示圓環(huán)邊緣的彈性變形范圍，根據(jù)實(shí)際情況通?？扇?．5甚至更小，即認(rèn)為點(diǎn)到圓心的距離誤差可以控制在0．5 m，這樣距離地面小于0．5 m的點(diǎn)可認(rèn)為是地面微小的起伏造成，被判斷為是地面細(xì)節(jié)予以保留。

彈性變形參數(shù)Δ固定，當(dāng)圓環(huán)半徑r較小時(shí)，可以滾過較小的凹凸起伏，保留住更多的地表細(xì)節(jié)，半徑較大時(shí)，則濾波結(jié)果比較平坦，較小的凹凸起伏會(huì)被濾除；圓環(huán)半徑r固定，當(dāng)彈性變形參數(shù)Δ較大時(shí)，可以保留更多細(xì)節(jié)，較小的Δ會(huì)濾掉更多細(xì)節(jié)。采用上述方法，某地形的掃描線濾波結(jié)果如圖6所示。

圖6 掃描線濾波結(jié)果Fig．6 Result of scanning line filtering

2．4 端點(diǎn)處理

如果地形的邊緣剛好是某建筑物或植被的最高點(diǎn)，滾圓法濾波后會(huì)出現(xiàn)端點(diǎn)上翹，即把端點(diǎn)處的地物點(diǎn)誤判為地面點(diǎn)而犯Ⅱ型錯(cuò)誤，如圖7（a）的左端點(diǎn)。此情況需要在端點(diǎn)附近取連續(xù)地面點(diǎn)擬合光滑曲線對(duì)端點(diǎn)進(jìn)行插值，因?yàn)榈匦尉植科鸱淮笸ǔ？刹捎枚味囗?xiàng)式插值，插值后如圖7（b）所示，端點(diǎn)回落到地面并避免了Ⅱ型錯(cuò)誤。

圖7 切面線濾波結(jié)果Fig．7 Result of section line filtering

3 B樣條曲面逼近

通過掃描線濾波，得到光滑的地形掃描線，并且地形掃描線之間間距相等。為了采用均勻有理B樣條曲面擬合地形表面，首先要在每條掃描線上等間隔采樣獲取型值點(diǎn)陣列，型值點(diǎn)是指位于地形表面上的點(diǎn)，即擬合后B樣條曲面上的點(diǎn)；其次，為了對(duì)原始點(diǎn)云中每個(gè)點(diǎn)進(jìn)行過濾，需要通過型值點(diǎn)陣列反算出擬合B樣條曲面的控制點(diǎn)陣列，控制點(diǎn)直接決定了擬合后B樣條曲面的形狀。最后，把原始點(diǎn)云中每個(gè)點(diǎn)的位置坐標(biāo)（x，y）代入B樣條參數(shù)曲面，求得點(diǎn)在樣條曲面上的投影高程，并與實(shí)際高程比較，根據(jù)給定的閾值判斷并區(qū)分地面點(diǎn)和地物點(diǎn)。

假設(shè)在等間隔的掃描線上均勻采樣得到型值點(diǎn)陣列P＝｛pij，i＝1，2，…，m，j＝1，2，…，n｝。即有m條等間隔掃描線，在每一條掃描線上均勻采樣n個(gè)點(diǎn)。由型值點(diǎn)陣列反算控制點(diǎn)陣列，可以分別在與地面平行的兩個(gè)方向上先后進(jìn)行兩次一維反算，通常要兩次求解式（11）所示的三對(duì)角線性方程組

式中，Pi表示型值點(diǎn)；Vi表示控制點(diǎn)；ai、bi和ci為線性方程組的系數(shù)，通過式（12）計(jì)算

式中，k表示曲線或曲面的階次。B樣條曲面的相關(guān)計(jì)算細(xì)節(jié)可以參考文獻(xiàn)［24］。

由于地形曲面不封閉，對(duì)于非周期曲線，控制點(diǎn)數(shù)比方程數(shù)多2，可補(bǔ)充切失邊界條件如式（14）即兩端點(diǎn)的切矢量的值給定分別為T0和Tn，可以由兩端型值點(diǎn)計(jì)算。

反算得到均勻有理B樣條曲面的控制點(diǎn)陣列，可以計(jì)算點(diǎn)云中任一點(diǎn)在樣條曲面上的投影，采取二維德布爾快速算法［24］計(jì)算該點(diǎn)在B樣條曲面上的投影點(diǎn)并計(jì)算其高程hij，比較投影點(diǎn)高程hij與點(diǎn)的實(shí)際高程zij，如果小于規(guī)定的閾值Δ，即滿足｜hij－zij｜＜Δ，該點(diǎn)為地面點(diǎn)予以保留，否則為地物點(diǎn)應(yīng)該濾除。閾值大小體現(xiàn)了地面允許的凹凸起伏程度，通?？梢匀?．5 m。閾值小，濾波后地形較光滑，過小則容易犯Ⅰ型錯(cuò)誤；閾值大，允許地面的起伏大，保留地面細(xì)節(jié)多，過大則容易犯Ⅱ型錯(cuò)誤。

4 試驗(yàn)與結(jié)果分析

4．1 試驗(yàn)數(shù)據(jù)

為了驗(yàn)證本文算法的有效性，本文采取ISPRS公布的數(shù)據(jù)進(jìn)行濾波試驗(yàn)，其中3個(gè)具有典型地貌特征的地形見表1。

表1 地形特征Tab．1 Terrain feature

4．2 試驗(yàn)結(jié)果定性分析

試驗(yàn)中采用的計(jì)算機(jī)配置為CPU 2．0 GHz，內(nèi)存2GB。3個(gè)典型地形濾波結(jié)果如圖8所示。圖8（a）、圖8（c）、圖8（e）為濾波前地形，圖8（b）、圖8（d）、圖8（f）分別為其對(duì)應(yīng)的濾波后地形，圖8（a）和圖8（c）是城市地形，圖8（e）為林區(qū)地形。對(duì)圖8（a）、圖8（c）、圖8（e）3個(gè)地形濾波的時(shí)間消耗分別為17．3、3．2和13．2 s，可見具有較高的時(shí)間效率。

圖8 濾波結(jié)果Fig．8 Results of filtering

從圖8中對(duì)不同地形的濾波結(jié)果看，疏密程度不同的地表植被和各種形狀的建筑物基本被濾除，同時(shí)保留了地形特征和地面細(xì)節(jié)。

4．3 試驗(yàn)結(jié)果定量分析

為便于對(duì)濾波算法進(jìn)行定量分析，文獻(xiàn)［1］提出了3個(gè)指標(biāo)：Ⅰ型誤差、Ⅱ型誤差和總誤差，計(jì)算方法為

Ⅰ型誤差＝誤判為地物的地面點(diǎn)數(shù)／地面點(diǎn)數(shù)Ⅱ型誤差＝誤判為地面的地物點(diǎn)數(shù)／地物點(diǎn)數(shù)總誤差＝總的誤判點(diǎn)數(shù)／總點(diǎn)數(shù)

采用ISPRS標(biāo)準(zhǔn)地形對(duì)算法進(jìn)行測試，限于篇幅本文僅列出3個(gè)典型地形測試結(jié)果，如圖9所示。

圖9 濾波結(jié)果Fig．9 Results of filtering

圖9（a）為覆蓋植被和建筑物的陡坡，圖9（c）為密集的建筑物地形，圖9（e）為變化頻繁的長地物地形。3個(gè)地形的濾波誤差如下表2—4所示。

表2 地形1濾波測試結(jié)果Tab．2 Filtering test result of terrain 1 （%）

表3 地形2濾波測試結(jié)果Tab．3 Filtering test result of terrain 2 （%）

表4 地形3濾波測試結(jié)果Tab．4 Filtering test result of terrain 3 （%）

試驗(yàn)表明，對(duì)于不同特征的地形，本文算法通用性好并且濾波誤差更小，從表2—4可以看出與傳統(tǒng)方法相比濾波精度提高了1～5倍。此外傳統(tǒng)方法需要迭代，例如Pfeifer提出的層次法、Axelsson提出的三角網(wǎng)漸進(jìn)加密法、Brovelli方法以及可變窗口的數(shù)學(xué)形態(tài)學(xué)方法等。迭代的方法最壞情況下其時(shí)間復(fù)雜度為O（n2），導(dǎo)致算法效率不高。本文算法不僅避免了迭代，而且通過獲取掃描線轉(zhuǎn)化為一維計(jì)算，其算法時(shí)間復(fù)雜度為O（n）。

5 結(jié) 論

本文提出了一種基于可變半徑圓環(huán)和B樣條擬合的地形點(diǎn)云濾波方法，該方法首先從地形表面面向地心的一側(cè)以可變半徑的滾圓法對(duì)等間隔的地形掃描線濾波，獲得光滑的地形表面采樣陣列，然后針對(duì)采樣陣列，采用均勻有理B樣條曲面進(jìn)行擬合，計(jì)算原始地形點(diǎn)云中每個(gè)點(diǎn)在擬合曲面上的投影點(diǎn)高程，并與實(shí)際高程比較，如果在給定閾值范圍內(nèi)，則被認(rèn)為是地面點(diǎn)并保留，否則被認(rèn)為是地物點(diǎn)并濾除。最后采取ISPRS標(biāo)準(zhǔn)數(shù)據(jù)進(jìn)行試驗(yàn)，驗(yàn)證了本文算法的有效性和通用性，同時(shí)也驗(yàn)證了本文算法的執(zhí)行速度和濾波精度。

［1］ SITHOLE G，VOSSELMAN G．Experimental Comparison of Filter Algorithms for Bare－Earth Extraction from Airborne Laser Scanning Point Clouds［J］．ISPRS Journal of Photogrammetry and Remote Sensing，2004，59（1－2）：85－101．

［2］ VOSSELMAN G．Slope Based Filtering of Laser Altimetry Data［J］．International Archives of Photogrammetry and Remote Sensing，2000，33（PartB 3／2）：935－942．

［3］ AXELSSON P．DEM Generation from Laser Scanner Data Using Adaptive TIN Models［J］．International Archives of the Photogrammetry and Remote Sensing，1999，33（B4／1）：110－117．

［4］ K ANG Xiaochen，LIU Jiping，LIN Xiangguo．Parallel Filter of Progressive TIN Densification for Airborne Li－DAR Point Cloud Using Multi－core CPU ［J］． Acta Geodaetica et Cartographica Sinica，2013，42（3）：331－336．（亢曉琛，劉紀(jì)平，林祥國．多核處理器的機(jī)載激光雷達(dá)點(diǎn)云并行三角網(wǎng)漸進(jìn)加密濾波方法［J］．測繪學(xué)報(bào)，2013，42（3）：331－336．）

［5］ LEE H S，YOUNAN N H．DTM Extraction of LiDAR Returns via Adaptive Processing［J］．IEEE Transactions on Geoscience and Remote Sensing，2003，41（9）：2063－2069．

［6］ WU Congcong，LU Xiaoping，LI Guoli，et al．Research on Filtering Algorithm for LiDAR Data Based on TIN［J］．Bulletin of Surveying and Mapping，2013（3）：32－35．（吳叢叢，盧小平，李國利，等．基于TIN的LiDAR數(shù)據(jù)濾波算法研究［J］．測繪通報(bào)，2013（3）：32－35．）

［7］ SUI Lichun，ZH A NG Yibin，ZH A NG Shuo，et al．Filtering of Airborne LiDAR Point Cloud Data Based on Progressive TIN［J］．Geomatics and Information Science of Wuhan University，2011，36（10）：1159－1163．（隋立春，張熠斌，張碩，等．基于漸進(jìn)三角網(wǎng)的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（10）：1159－1163．）

［8］ KRAUS K，PFEIFER N．Determination of Terrain Models in Wooded Areas with Airborne Laser Scanner Data［J］．ISPRS Journal of Photogrammetry and Remote Sensing，1998，53（4）：193－203．

［9］ ELMQVIST M，Jungert E，Persson L?，et al．Terrain Modelling and Analysis Using Laser Scanner Data［J］．International Archives of the Photogrammetry and Remote Sensing，2001，34（3／W4）：219－224．

［10］ BROVELLI M A，CANNATA M．Digital Terrain Model Reconstruction in Urban Areas from Airborne Laser Scanning Data：the Method and the Example of the Town of Pavia（Northern Italy）［J］．International Archives of the Photogrammetry and Remote Sensing，2002．34（2）：43－48．

［11］ SH ANG Dashuai，LI Hengming，ZH AO Xi．A Filtering Method Based on Surface Constraint for Point Cloud Data［J］．Radar Science and Technology，2013，11（4）：413－418．（尚大帥，黎恒明，趙羲．一種基于曲面約束的點(diǎn)云數(shù)據(jù)濾波方法［J］．雷達(dá)科學(xué)與技術(shù)，2013，11（4）：413－418．）

［12］ KILIAN J，H AALA N，ENGLICH M．Capture and Evaluation of Airborne Laser Scanner Data ［J］．International Archives of the Photogrammetry and Remote Sensing，1996，31（B3）：383－88．

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

［14］ SUI Lichun，YANG Yun．Filtering of Airborne LiDAR Point Cloud Data Based on car（p，q） Model and Mathematical Morphology ［J］． Acta Geodaetica et Cartographica Sinica，2012，41（2）：219－224．（隋立春，楊耘．基于car（p，q）模型和數(shù)學(xué)形態(tài)學(xué)理論的LiDAR點(diǎn)云數(shù)據(jù)濾波［J］．測繪學(xué)報(bào)，2012，41（2）：219－224．）

［15］ ZHANG Keqi，CHEN Shuching，WHITMAN D，et al．A Progressive Morphological Filter for Removing Nonground Measurements from Airborne LiDAR Data［J］．IEEE Transactions on Geoscience and Remote Sensing，2003，41（4）：872－882．

［16］ CHEN Qi，GONG Peng，BALDOCCHI D，et al．Filtering Airborne Laser Scanning Data with Morphological Methods［J］．Photogrammetric Engineering and Remote Sensing，2007，73（2）：175－185．

［17］ SU N Meiling，LI Yongshu，CH EN Qiang，et al．A Progressive Morphological Filtering Method for Airborne LiDAR Point Cloud Based on Scan Line［J］．Opto－Electronic Engineering，2013，40（11）：71－75．（孫美玲，李永樹，陳強(qiáng)，等．基于掃描線的漸進(jìn)式形態(tài)學(xué)機(jī)載LiDAR點(diǎn)云濾波［J］．光電工程，2013，40（11）：71－75．）

［18］ SHEN Jing，LIU Jiping，LIN Xiangguo．Airborne LiDAR Data Filtering by Morphological Reconstruction Method［J］．Geomatics and Information Science of Wuhan University，2011，36（2）：167－170，175．（沈晶，劉紀(jì)平，林祥國．用形態(tài)學(xué)重建方法進(jìn)行機(jī)載LiDAR數(shù)據(jù)濾波［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2011，36（2）：167－170，175．）

［19］ LI Yong，WU Huayi．Filtering Airborne LiDAR Data Based on Morphological Gradient［J］．Journal of Remote Sensing，2008，12（4）：633－639．（李勇，吳華意．基于形態(tài)學(xué)梯度的機(jī)載激光掃描數(shù)據(jù)濾波方法［J］．遙感學(xué)報(bào)，2008，12（4）：633－639．）

［20］ TANG Feifei，LIU Jingnan，ZHANG Xiaohong，et al．A Voxel－based Filtering Algorithm For DTM Data Extraction in Forest Areas［J］．Journal of Beijing Forestry University，2009，31（1）：55－59．（唐菲菲，劉經(jīng)南，張小紅，等．基于體素的森林地區(qū)機(jī)載LiDAR數(shù)據(jù)DTM提取［J］．北京林業(yè)大學(xué)學(xué)報(bào)，2009，31（1）：55－59．）

［21］ ZH A NG Xiaohong．Airborne LiDAR Points Cloud Classification Based on Difference of Twice Return Pulse Heights［J］．Science of Surveying and Mapping，2006，31（4）：48－50．（張小紅．利用機(jī)載LiIDAR雙次回波高程之差分類激光腳點(diǎn)［J］．測繪科學(xué)，2006，31（4）：48－50．）

［22］ WAN Jianhua，HUANG Ronggang，ZHOU Hang，et al．A Secondary Filter Method of LiDAR Point Cloud Based on Curvature Statistics［J］．Journal of China University of Petroleum，2013，37（1）：56－60．（萬劍華，黃榮剛，周行，等．基于曲率統(tǒng)計(jì)的LiDAR點(diǎn)云二次濾波方法［J］．中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，37（1）：56－60．）

［23］ HUANG Xianfeng，LI Hui，WANG Xiao，et al．Filter Algorithms of Airborne LiDAR Data：Review and Prospects［J］．Acta Geodaetica et Cartographica Sinica，2009，38（5）：466－469．（黃先鋒，李卉，王瀟，等．機(jī)載LiDAR數(shù)據(jù)濾波方法評(píng)述［J］．測繪學(xué)報(bào)，2009，38（5）：466－469．）

［24］ SHI Fazhong．CAGD ＆ NURBS［M］．2nd ed．Beijing：Higher Education Press，2013．（施法中．計(jì)算機(jī)輔助幾何設(shè)計(jì)與非均勻有理B樣條［M］．第2版．北京：高等教育出版社，2013．）