萬年鋒 袁啟倫 邸國輝 周勝潔

摘要：

為了解決傳統(tǒng)的不規(guī)則三角網(wǎng)（Triangular irregular network，TIN）算法在陡崖反坡地形濾波時，由于陡崖反坡會產(chǎn)生多個投影交點，坡度出現(xiàn)負值，陡崖反坡區(qū)域巖壁點往往被錯誤過濾的問題，提出一種剛性坐標(biāo)變換與TIN結(jié)合的機載激光雷達（Light detection and ranging，LiDAR）點云濾波算法。采用剛性坐標(biāo)變換的方法重新構(gòu)建陡崖面的直觀形態(tài)，將陡崖反坡等特異地形常規(guī)化。以湖北省恩施市姚家平水利樞紐工程壩址區(qū)域的陡崖反坡數(shù)據(jù)進行濾波分析，構(gòu)建了側(cè)向數(shù)字高程模型（DEM）表達陡崖地形，實現(xiàn)了側(cè)向DEM在壩體填筑量計算上的應(yīng)用。結(jié)果表明：該算法能夠在陡崖反坡區(qū)域大幅降低Ⅰ類誤差，并將Ⅱ類誤差控制在一定范圍內(nèi)；應(yīng)用該技術(shù)可為大壩設(shè)計提供直觀準(zhǔn)確的DEM，為三維設(shè)計提供了新思路。

關(guān)鍵詞：

陡崖反坡地形； 點云濾波； 機載LiDAR； 不規(guī)則三角網(wǎng)； 側(cè)向DEM； 大壩設(shè)計

中圖法分類號：P237

文獻標(biāo)志碼：A

DOI：10.15974/j.cnki.slsdkb.2023.05.007

文章編號：1006-0081（2023）05-0045-05

0 引 言

機載激光雷達（Light detection and ranging，LiDAR）技術(shù)結(jié)合了全球定位系統(tǒng)、激光掃描儀與慣導(dǎo)系統(tǒng)，能夠快速高效地獲取高精度地面點云數(shù)據(jù)［1］，具有效率高、受影響小、精度高、風(fēng)險小等優(yōu)勢［2-4］。由于原始點云是由不同地形特征（地面、建筑物等）返回的大量離散點組合而成，并無語義信息或拓撲關(guān)聯(lián)，因此在應(yīng)用之前，必須將地面和非地面點分離，這個過程稱為點云濾波。Lindenberger［5］較早應(yīng)用數(shù)學(xué)形態(tài)學(xué)對剖面式測高數(shù)據(jù)執(zhí)行開運算，然后使用自回歸過程調(diào)整運算結(jié)果，實現(xiàn)對原始數(shù)據(jù)濾波，該算法原理簡單、實現(xiàn)效率高，但缺乏適應(yīng)性，容易削平地形。Axelsson［6］提出了三角網(wǎng)加密濾波算法，采用不規(guī)則三角形擬合地形，不斷將滿足地形特征的點加入三角網(wǎng)，通過三角網(wǎng)加密的形式還原地形信息，該算法適應(yīng)性較強，已被應(yīng)用于濾波處理軟件TerraSolid模塊中，但易丟失地面點且低矮植被不能有效剔除。Kraus等［7］提出了最小二乘內(nèi)插濾波算法，根據(jù)待判點與擬合曲面之間的高程差確定相應(yīng)的權(quán)重，然后依據(jù)權(quán)重值對擬合曲面進行調(diào)整，反復(fù)以上操作，直到完成所有數(shù)據(jù)，該算法濾波精度較高，但需要多次迭代，運行速度慢。羅伊萍等［8］采用多尺度的方式有效改善了固定窗口大小，提升了數(shù)學(xué)形態(tài)學(xué)算法的適應(yīng)性，但組織形式有損原始數(shù)據(jù)的精度。徐國杰等［9］采用設(shè)定點到三角面的最大距離和點與三角面頂點的連線形成的最大角的方式對不規(guī)則三角網(wǎng)進行改進，并對參數(shù)動態(tài)調(diào)整，多次迭代完成濾波，但運行效率較差。張皓等［10］使用虛擬網(wǎng)格組織數(shù)據(jù)，提出了一種改進坡度的濾波算法，可以有效避免降低數(shù)據(jù)精度，但閾值的自適應(yīng)性較差。上述各種主流的LiDAR濾波與分類算法研究都有其各自特點及局限性。其中，三角網(wǎng)濾波算法在山區(qū)等坡度變化大的地區(qū)適應(yīng)性最佳，總誤差值相對較小，濾波時間也更短［11］，但是對于山區(qū)里的一些特殊地形，比如水利工程中常見的陡崖及反崖地形，運用三角網(wǎng)濾波算法時，會出現(xiàn)一些陡崖崖壁及反崖反向坡度的數(shù)據(jù)損失，無法真實表達其地理地貌［12］。本文首次提出一種剛性坐標(biāo)變換與不規(guī)則三角網(wǎng)（triangular irregular network，TIN）相結(jié)合的機載LiDAR點云濾波算法，能在陡崖及反坡等特殊地形高效準(zhǔn)確進行濾波的方法，并結(jié)合湖北省恩施市姚家平水利樞紐工程，對無人機載LiDAR在陡崖反坡地形測繪及大壩設(shè)計中的應(yīng)用進行實際分析。

1 陡崖反坡濾波方法

1.1 陡崖反坡

目前LiDAR點云濾波軟件均未考慮陡崖反坡存在的多個投影交點（如圖1所示的 O，P，Q ）、坡度出現(xiàn)負值的情況，在點云濾波提取地面點時，陡崖區(qū)域巖壁點（ P和Q ）往往被忽略。本文提出了對陡崖反坡點云先進行變換處理再進行濾波的方法。

1.2 陡崖反坡點云變換處理

在掃描獲取陡崖點云數(shù)據(jù)后，采用高精度定位定姿系統(tǒng)（Position and orientation system，POS）數(shù)據(jù)解算和激光點云數(shù)據(jù)解算法對獲取的數(shù)據(jù)進行解析，得到工程坐標(biāo)系（ XYH ）0的點云，并構(gòu)建陡崖區(qū)域，剪裁后，得到陡崖區(qū)域點云文件LAS0。工程坐標(biāo)系（ XYH ）0為項目工程指定的坐標(biāo)系，該坐標(biāo)系包 括X軸、Y軸和H軸 。

基于已有數(shù)據(jù)繪制出整個陡崖的上邊緣線，并對上邊緣線進行直線擬合，得到線段 AB，作法線垂直于AB。AB的 法線方向如圖2所示。對點云LAS0進行的剛性坐標(biāo)變換，包括第一次剛性變換和第二次剛性變換。

（1） 第一次剛性變換。對點云LAS0進行第一次剛性變換，將其縱坐標(biāo)軸旋轉(zhuǎn) θ角至線段AB 的法線方向，與法線方向重合（圖2），得到坐標(biāo)系（ XYH ）1，其中坐標(biāo)按式（1）換算：

x ′ y′ h′= cos θ－ sin θ0? sin θ cos θ0 001x y h（1）

式中： x y h為初始點云坐標(biāo)；x′ y′ h′為第一次剛性變換后點云坐標(biāo)；θ為 旋轉(zhuǎn)角。

（2） 第二次剛性變換。將第一次剛性變換后的高程坐標(biāo)作為新的 X軸坐標(biāo)，第一次剛性變換后的X軸坐標(biāo)轉(zhuǎn)換為新的Y軸坐標(biāo)，第一次剛性變換后的Y軸坐標(biāo)轉(zhuǎn)換為新的高程坐標(biāo)，得到新坐標(biāo)系下陡崖面的初始點云LAS1。以p為原始點云上的任一點，p′為第一次剛性變換后p點的對應(yīng)點，p″為第二次剛性變換后p點的對應(yīng)點。p點坐標(biāo)為（x，y，h），p′坐標(biāo)為（x′，y′，h′），p″坐標(biāo)為（x″，y″，h″）。p繞原點旋轉(zhuǎn)θ，得到點p′，相應(yīng)的x，y坐標(biāo)也旋轉(zhuǎn)θ為x′，y′，高程h′等于h。再對p′進行第二次剛性變換：高程h′坐標(biāo)值轉(zhuǎn)換為第二次剛性變換后的X軸，即X″；X′坐標(biāo)值轉(zhuǎn)換為第二次剛性變換后的Y軸，即Y″；Y′坐標(biāo)值轉(zhuǎn)換為第二次剛性變換后的H軸，即H″，得到坐標(biāo)系（XYH） 2的點云，即陡崖面的初始點云。坐標(biāo)換算如式（2）所示：

x″ y″ h″=001 100 010x′ y′ h′（2）

第二次剛性變換后， x″的數(shù)值等于原高程h，而x″ 作為橫坐標(biāo)，投影方向與法線方向一致，消除了陡崖反坡形態(tài)，顯著減小了陡崖面坡度，得到陡崖面的初始點云LAS1，可直觀反映陡崖的形態(tài)。此外，由于消除了反坡形態(tài)，深度坐標(biāo)面的坡度顯著減小，實現(xiàn)了陡崖區(qū)域巖壁地表信息的高效提取。至此，LAS1的反坡被消除。

1.3 濾 波

利用點云處理軟件TerraSolid，基于不規(guī)則三角網(wǎng)TIN加密濾波算法，先通過一些較低種子點生成一個稀疏TIN，然后考察每個點與TIN的距離和坡度，并逐層迭代加密，進行地面點和非地面點（植被）分類，由LAS1提取出地面點點云LAS2。

2 實例分析

2.1 研究區(qū)域概況

研究區(qū)位于湖北省恩施市姚家平水利樞紐，試驗航攝區(qū)域面積約1 km2，測區(qū)平均海拔700 m，最大相對高差為800 m。測區(qū)沿河谷方向分布多處陡崖及反崖。

2.2 無人機LiDAR航飛

航飛采用華測BB4大黃蜂四旋翼無人機+AS-900HL多平臺激光雷達系統(tǒng)，精度高、測程長，配備高精度光纖慣導(dǎo)系統(tǒng)，滿足此次LiDAR掃描任務(wù)。

由于測區(qū)屬于峽谷地區(qū)，高差大，故整個測區(qū)采用仿地飛行。仿地飛行需要的測區(qū)初始數(shù)字地表模型（Digital surface model，DSM）來源于前期0.2 m分辨率的傾斜攝影成果。為掃描陡崖和反坡區(qū)域的地形，在平行峽谷方向布設(shè)了4條航線。航線高度120 m，航向重疊率60%，旁向重疊均值35%，點云密度不少于8個/m2。測區(qū)初始DSM及仿地飛行航線如圖3所示。

2.3 陡崖及反崖濾波

沿北岸峽谷方向的谷底有1條200 m長的陡崖反坡，其原始剖面點云見圖4。

（1） 基于1.2節(jié)的剛性變換，得到新坐標(biāo)系下陡崖面的初始點云LAS1，如圖5所示。可明顯觀察到，變換處理消除了陡崖反坡形態(tài)，顯著減小了陡崖面坡度。

（2） 利用點云處理軟件TerraSolid，基于不規(guī)則三角網(wǎng)TIN加密濾波算法，進行地面點和非地面點（植被）分類，由LAS1提取出地面點點云LAS2。其中，反坡設(shè)置參數(shù)：坡度80°，建筑物長度10 m，反復(fù)角設(shè)置為12°，提取地面點。

至此，得到了經(jīng)過空間轉(zhuǎn)換后的濾波點云。若需要得到原坐標(biāo)系下濾波點云，經(jīng)過反向變換即可。此外，由于陡崖（特別是反崖）的特殊性，在構(gòu)建DEM時仍會出現(xiàn)因為DEM數(shù)據(jù)水平面上某一點只能對應(yīng)唯一的高程值而無法生成反坡處DEM的問題，導(dǎo)致陡崖反坡的真實地形無法表達。為此，提出構(gòu)建側(cè)向DEM的方法表達陡崖地形，更能滿足實際工程設(shè)計的需要。

（3） 側(cè)向DEM是利用Globalmapper軟件，導(dǎo)入陡崖面地面點的點云LAS2生成的DEM，如圖6所示。因為點云LAS2是經(jīng)過變換旋轉(zhuǎn)后的點云，數(shù)學(xué)關(guān)系如1.2節(jié)的點云變換原理所述，故稱為側(cè)向DEM。這種經(jīng)過變換旋轉(zhuǎn)后的DEM可以更好地表達反坡處地形，有利于后期設(shè)計。

2.4 結(jié)果分析

根據(jù)2003年國際攝影測量和遙感學(xué)會（ISPRS）發(fā)布的誤差評價標(biāo)準(zhǔn)，對本次實驗濾波結(jié)果進行定量評價［13］。此種誤差評價方法根據(jù)激光點的誤分率，大致可以分為3類：① Ⅰ類誤差，即將地面點分類為非地面點的誤差（拒真誤差，omission error）；② Ⅱ類誤差，即將非地面點分類為地面點的誤差（納偽誤差，commission error）；③ Ⅲ類誤差，即總誤差。前兩類反映了算法的適應(yīng)性，總誤差反映了算法的可行性。為了最大限度得到地形表面數(shù)據(jù)，應(yīng)盡可能避免Ⅱ類誤差［14］。

為便于對比分析，同時采用傳統(tǒng)TIN算法及本文算法對陡崖反坡點云數(shù)據(jù)進行濾波實驗和精度評價，結(jié)果見表1。與傳統(tǒng)TIN算法相比，本文算法的Ⅰ類誤差降低幅度最大，降低34.02％；Ⅱ類誤差值大致不變；總誤差降低16.65％。綜上所述，本文算法在Ⅱ類誤差無較大變化的同時，顯著降低Ⅰ類誤差，驗證了該算法具有較好的普適性與有效性，能夠獲得更接近真實地形的高精度數(shù)字高程模型。

2.5 重力壩填筑量計算

結(jié)合恩施市姚家平水利樞紐工程的重力壩方案，采用側(cè)向DEM計算壩體填筑量［15］。取LAS0的縱軸至壩體軸線方向，將壩體順時針旋轉(zhuǎn)90°，將壩體橫斷面導(dǎo)入側(cè)向DEM，給定基點高程（左或右壩段長度），進行填筑量計算。

壩高97 m，內(nèi)坡比1∶0.25，外坡比1∶0.15，底部寬度53.3 m，頂部寬度14.5 m，取基點高程425 m（相應(yīng)左岸壩段長60 m），左岸壩段TIN地面點有8 661個，點密度為2.6個/m2。采用DEM棱柱法計算方量，得填筑量為210 489.68 m3，見圖7。同理，右岸壩段也可通過這種方法計算，進而得到整個壩段的填筑量。

通常斷面法土方量的計算采用由兩相鄰橫斷面的平均面積乘以兩相鄰橫斷面間距的方法（斷面法土方計算平均公式），即由側(cè)向DEM截取11個斷面計算斷面面積，計算所得填筑量為209 331.98 m3，與DEM棱柱法對比，相對誤差0.55%。結(jié)果表明：DEM棱柱法計算方量的誤差遠小于SL 197-97《水電水利工程施工測量規(guī)范》的要求。

DEM棱柱法充分利用了陡崖面地面點的點云數(shù)據(jù)，計算方法嚴(yán)密，比斷面法計算精度高，為三維設(shè)計提供了新方法。

3 結(jié) 論

（1） 本文在TIN算法基礎(chǔ)上，提出了通過對點云進行空間變換，改變陡崖及反坡的直觀形態(tài)的方法，改善了濾除效果和地形起伏區(qū)域的適應(yīng)性，顯著提高了濾波精度。與傳統(tǒng)TIN算法相比，該算法能有效過濾出山區(qū)陡崖及反坡的真實地面點。

（2） 利用恩施市姚家平水利樞紐工程區(qū)域陡崖反坡機載LiDAR點云數(shù)據(jù)進行實例分析，其試驗結(jié)果定量證明了剛性坐標(biāo)變換與TIN結(jié)合的機載LiDAR點云濾波算法的精度及實用性。

（3） 在大壩設(shè)計方面，應(yīng)用該方法生成側(cè)向DEM進行重力壩填筑量計算，結(jié)果表明：DEM棱柱法計算方量的誤差小于規(guī)范要求，可用于大壩三維設(shè)計。

參考文獻：

［1］ 易志朝.機載LiDAR技術(shù)在水利工程大比例尺地形圖測繪中的應(yīng)用［J］.陜西水利，2022（5）：181-183.

［2］ 談?wù)Z，趙云毅，富博.復(fù)雜山區(qū)河道地形機載Lidar數(shù)據(jù)質(zhì)量控制［J］.水利水電快報，2021，42（10）：22-26.

［3］ 劉世振，鄧建華，馮國正，等.機載LiDAR在山區(qū)型河道地形測繪中的適用性研究［J］.人民長江，2021，52（1）：108-113.

［4］ 吳昊，張曉萌.機載LiDAR在長江中游河道測繪中的應(yīng)用［J］.人民長江，2022，53（9）：109-114.

［5］ LINDENBERGER J.Laser-Profilmessungen zur topographischen Gelndeaufnahme［D］.Stuttgart：Stuttgart University，1993.

［6］ AXELSSON P.DEM generation from laser scanner data using adaptive TIN models［J］.International Archives of Photogrammetry & Remote Sensing，2000，23（B4）：110-117.

［7］ KRAUS K，PFEIFER N.Determination of terrain models in wooded areas with airborne laser scanner data［J］.ISPRS Journal of Photogrammetry & Remote Sensing，1998，53（4）：193-203.

［8］ 羅伊萍，姜挺，龔志輝，等.基于自適應(yīng)和多尺度數(shù)學(xué)形態(tài)學(xué)的點云數(shù)據(jù)濾波方法 ［J］.測繪科學(xué)技術(shù)學(xué)報，2009，26（6）：426-429.

［9］ 徐國杰，胡文濤.一種改進的 LiDAR 點云 TIN 迭代濾波算法［J］.測繪地理信息，2010，35（1）：33-35.

［10］ 張皓，賈新梅，張永生，等.基于虛擬網(wǎng)格與改進坡度濾波算法的機載LIDAR數(shù)據(jù)濾波［J］.測繪科學(xué)技術(shù)學(xué)報，2009，26（3）：224-227，231.

［11］ 鄒正，鄒進貴，胡海洋.不同機載LiDAR點云濾波算法對比分析［J］.測繪地理信息，2021，46（5）：52-56.

［12］ 朱依民，田林亞，畢繼鑫，等.基于PTD和改進曲面擬合的高山區(qū)水電工程機載激光雷達點云濾波方法［J］.水利水電科技進展，2021，41（1）：35-40.

［13］ 吉雨田，張春亢，尹耀.機載LiDAR點云自適應(yīng)濾波算法［J］.測繪科學(xué)技術(shù)學(xué)報，2021，38（2）：142-147.

［14］ 馬超，蔣好忱，秦先鋒，等.LiDAR-DP點云自動分類及精度評定［J］.測繪與空間地理信息，2022，45（4）：240-242，246.

［15］ 孟永東，陳俊華.約束TIN在土石方工程三維輔助設(shè)計中的應(yīng)用［J］.計算機應(yīng)用，2010，30（增1）285-288.

（編輯：江 燾，高小雲(yún)）

Abstract：

In order to solve the problem that the steep cliff reverse slope will produce multiple projection intersections，the slope will have negative values，and the rock wall points in the steep cliff reverse slope area are often wrongly filtered when the traditional triangular irregular network （TIN） algorithm is used to filter the terrain of the steep cliff reverse slope，a light detection and ranging（LiDAR） point cloud filtering algorithm combined with rigid coordinate transformation and TIN was proposed.The intuitive shape of the steep cliff surface was reconstructed by rigid coordinate transformation，which could normalize the specific terrain such as the steep cliff reverse slope to facilitate further filtering.Lateral Digital Elevation Model（DEM） was constructed to express the steep cliff terrain，the application of lateral DEM on dam filling amount calculation was realized.The results showed that this algorithm could greatly reduce the class Ⅰ error in the steep cliff reverse slope area and control the class Ⅱ error within a certain range.This method could provide intuitive and accurate DEM for dam design，which provided a new method for three-dimensional design.

Key words：

steep cliff reverse slope； point cloud filtering； UAV LiDAR； triangular irregular network； lateral DEM； dam design