基于無人機(jī)LiDAR數(shù)據(jù)的陡崖反坡濾波技術(shù)及應(yīng)用

萬(wàn)年鋒 袁啟倫 邸國(guó)輝 周勝潔

摘要：

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

關(guān)鍵詞：

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

中圖法分類號(hào)：P237

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

DOI：10.15974/j.cnki.slsdkb.2023.05.007

文章編號(hào)：1006-0081（2023）05-0045-05

0 引 言

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

1 陡崖反坡濾波方法

1.1 陡崖反坡

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

1.2 陡崖反坡點(diǎn)云變換處理

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

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

（1） 第一次剛性變換。對(duì)點(diǎn)云LAS0進(jìn)行第一次剛性變換，將其縱坐標(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為初始點(diǎn)云坐標(biāo)；x′ y′ h′為第一次剛性變換后點(diǎn)云坐標(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)系下陡崖面的初始點(diǎn)云LAS1。以p為原始點(diǎn)云上的任一點(diǎn)，p′為第一次剛性變換后p點(diǎn)的對(duì)應(yīng)點(diǎn)，p″為第二次剛性變換后p點(diǎn)的對(duì)應(yīng)點(diǎn)。p點(diǎn)坐標(biāo)為（x，y，h），p′坐標(biāo)為（x′，y′，h′），p″坐標(biāo)為（x″，y″，h″）。p繞原點(diǎn)旋轉(zhuǎn)θ，得到點(diǎn)p′，相應(yīng)的x，y坐標(biāo)也旋轉(zhuǎn)θ為x′，y′，高程h′等于h。再對(duì)p′進(jìn)行第二次剛性變換：高程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的點(diǎn)云，即陡崖面的初始點(diǎn)云。坐標(biāo)換算如式（2）所示：

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

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

1.3 濾 波

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

2 實(shí)例分析

2.1 研究區(qū)域概況

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

2.2 無人機(jī)LiDAR航飛

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

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

2.3 陡崖及反崖濾波

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

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

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

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

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

2.4 結(jié)果分析

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

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

2.5 重力壩填筑量計(jì)算

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

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

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

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

3 結(jié) 論

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

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

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

參考文獻(xiàn)：

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

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

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

［4］ 吳昊，張曉萌.機(jī)載LiDAR在長(zhǎng)江中游河道測(cè)繪中的應(yīng)用［J］.人民長(zhǎng)江，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é)的點(diǎn)云數(shù)據(jù)濾波方法 ［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2009，26（6）：426-429.

［9］ 徐國(guó)杰，胡文濤.一種改進(jìn)的 LiDAR 點(diǎn)云 TIN 迭代濾波算法［J］.測(cè)繪地理信息，2010，35（1）：33-35.

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

［11］ 鄒正，鄒進(jìn)貴，胡海洋.不同機(jī)載LiDAR點(diǎn)云濾波算法對(duì)比分析［J］.測(cè)繪地理信息，2021，46（5）：52-56.

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

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

［14］ 馬超，蔣好忱，秦先鋒，等.LiDAR-DP點(diǎn)云自動(dòng)分類及精度評(píng)定［J］.測(cè)繪與空間地理信息，2022，45（4）：240-242，246.

［15］ 孟永東，陳俊華.約束TIN在土石方工程三維輔助設(shè)計(jì)中的應(yīng)用［J］.計(jì)算機(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