王喆，劉鵬飛，王凱

(北京市勘察設(shè)計研究院有限公司，北京 100038)

1 引 言

三維激光掃描技術(shù)(3D Laser Scanning Technology)作為一項高新測繪技術(shù)興起于20世紀90年代。相比于傳統(tǒng)測量手段，該技術(shù)采用非接觸的方式連續(xù)、自動采集目標對象的表面屬性點信息，數(shù)據(jù)形式為三維點云[1]。目前，三維激光掃描技術(shù)已在多個領(lǐng)域，例如城市科學規(guī)劃、設(shè)計、管理中普遍應(yīng)用。

傳統(tǒng)的地質(zhì)研究中，通常使用實地勘測的方法獲取巖體中結(jié)構(gòu)面的信息，由于該方法依靠技術(shù)人員以羅盤、皮尺等工具進行實地測量，故受現(xiàn)場局限性較大。當受到場地限制(如高陡邊坡)，通常難以獲取全面的斷面信息。由于此類基礎(chǔ)斷面數(shù)據(jù)不齊全，會對進一步的巖體裂隙分布規(guī)律判斷產(chǎn)生干擾，最終影響巖體穩(wěn)定性分析評價工作。并且此方法存在工作效率低、精度低、成本高、外業(yè)工作量大和安全隱患大等問題[2]。

利用三維激光掃描技術(shù)可實現(xiàn)無接觸地直接采集目標巖體的點云數(shù)據(jù)，數(shù)據(jù)經(jīng)過處理后得到采集目標的三維點云模型。由于該模型真實地逼近于現(xiàn)場情況，技術(shù)人員可在后期處理軟件中對危巖體幾何參數(shù)進行量測和計算，得到巖體結(jié)構(gòu)面產(chǎn)狀，從而完成一系列地質(zhì)分析。故此方法在地質(zhì)勘察工作中，尤其是大高差山區(qū)危巖地區(qū)，可基本代替?zhèn)鹘y(tǒng)人工勘測方法，對地質(zhì)體信息提取工作具有極大的價值[3]。

2 PCL簡介

PCL(Point Cloud Library)是一個大型跨平臺(Windows、Linux、MacOSX、Android等)C++編程庫，作為一個開源項目，其最早是由斯坦福大學的Radu博士等人維護和開發(fā)[4]。PCL是基于Boost、FLANN、Eigen、VTK、CUDA、OpenNI等第三方庫集成的模塊化C++模板庫?；赑CL可建立高效的點云數(shù)據(jù)結(jié)構(gòu)，依此結(jié)構(gòu)處理點云大數(shù)據(jù)，可實現(xiàn)相關(guān)的通用算法，PCL均將其進行了封裝并提供調(diào)用接口，其中包括點云的獲取、分割、檢索、配準、濾波、可視化、特征提取等。圖1為PCL的架構(gòu)圖[5]。

圖1 PCL架構(gòu)圖

2.1 1PCD文件

PCL作為針對點云數(shù)據(jù)開發(fā)的專業(yè)庫，為了高效地處理點云大數(shù)據(jù)，此庫采用了C++作為底層編程語言，并以此發(fā)布了PCL獨有的點云文件數(shù)據(jù)格式：PCD格式。隨著三維激光掃描技術(shù)在各行業(yè)領(lǐng)域的廣泛應(yīng)用，技術(shù)設(shè)備不斷迭代更新，各大設(shè)備和軟件廠商紛紛推出了各自專屬的點云數(shù)據(jù)格式，其中較為主流的格式包括LAS、OBJ、PLY、X3D、STL等。但上述數(shù)據(jù)格式大都是根據(jù)本廠家特定的應(yīng)用場景，在舊式傳感器設(shè)備和算法基礎(chǔ)上創(chuàng)建的。相比以上數(shù)據(jù)格式，PCD格式的主要優(yōu)勢如下[6]：

(1)采用mmap、munmap作為基礎(chǔ)二進制數(shù)據(jù)類型，使其擁有高效的數(shù)據(jù)下載和存儲能力。

(2)可存儲多類C++基本數(shù)據(jù)類型，如int、float、double、char、short等，并針對無效的數(shù)據(jù)點設(shè)置特定的數(shù)據(jù)類型：NAN類型，以便于在PCL內(nèi)部存儲，從而提高PCL的點云處理能力。

(3)通過創(chuàng)建PCD文件格式，可以發(fā)揮C++強大的數(shù)據(jù)處理能力的優(yōu)勢，并且能最大程度上適應(yīng)PCL與PCL的內(nèi)部函數(shù)無縫銜接，從而讓基于PCL開發(fā)的應(yīng)用程序獲得最佳的數(shù)據(jù)處理性能。

2.2 PCD文件操作

基于C++強大的底層支撐，PCL內(nèi)部提供了豐富的點云文件操作模塊，以模塊提供的一系列接口開發(fā)程序，可以很輕松地實現(xiàn)針對點云的相關(guān)操作。以 I/O模塊為例，該模塊主要用來支持PCD文件的輸入輸出操作，共包含21個子類。例如PCD文件讀寫的接口就可以使用File-Reader類和File-Writer類分別定義。二者的繼承關(guān)系如圖2所示，其功能實現(xiàn)的關(guān)鍵代碼為：

pcl：：PointCloud：：Ptrcloud(newpcl：：PointCloud)；

pcl：：io：：loadPCDFile(“test_pcd.pcd”，*cloud)；

pcl：：io：：savePCDFileASCII(“test_pcd.pcd”，cloud)；

圖2 類File Reader和File-Writer的繼承關(guān)系

3 平面分割算法

利用三維激光掃描獲取的點云數(shù)據(jù)，具有連續(xù)、實時、數(shù)據(jù)量大，精度高等特點，經(jīng)過矯正、降噪等預(yù)處理操作后，點云數(shù)據(jù)仍是一個龐大的完整巖體[7]。故欲進行地質(zhì)體產(chǎn)狀要素的計算與統(tǒng)計，需先將每一個結(jié)構(gòu)面單獨提出。本文采用RANSAC算法作為結(jié)構(gòu)面提取的算法，RANSAC算法最早由Fischler和Bolles于1981年提出，該算法以一組包含異常數(shù)據(jù)的樣本數(shù)據(jù)集為依據(jù)，計算出數(shù)據(jù)的模型參數(shù)，從而得到有效樣本[8]。其基本思想是：每一個平面數(shù)據(jù)由該平面的“局內(nèi)點”組成，其數(shù)據(jù)的分布可滿足于特定的模型解釋參數(shù)；刨除“局內(nèi)點”，數(shù)據(jù)源中的其他點再分為“局外點”與“噪聲點”，其中“局外點”是不能適應(yīng)“局內(nèi)點”模型的數(shù)據(jù)，剩余數(shù)據(jù)即為“噪聲點”。該算法結(jié)果的合理性是有一定概率限制的，因此提高迭代次數(shù)才能將結(jié)果合理性的概率提高，即要保證在一定置信度下基本子集最小抽樣數(shù)N與至少取得一個良好抽樣子集的概率P滿足如下關(guān)系[9]：

P=1-(1-εk)n

(1)

式中，k為計算模型參數(shù)需要的最小數(shù)據(jù)量；ε為局內(nèi)點與數(shù)據(jù)集點數(shù)的比值；P一般取值為0.9～0.99。對上式兩邊取對數(shù)可得：

(2)

算法迭代的核心邏輯為：

(1)首先隨機選擇一個面為最佳擬合平面，比較當前迭代的點數(shù)與該最佳平面點數(shù)，如當前迭代點數(shù)更多，將當前點設(shè)置為最佳擬合平面；

(2)如當前點數(shù)量與最佳擬合平面點數(shù)量相同，但當前點滿足更小閾值時，同樣將當前點設(shè)置為最佳擬合平面。

總體上，RANSAC對噪聲點有一定的抑制作用，并能較好地分割多面物體，是一種有效的穩(wěn)健估計算法[9]。

具體算法流程如下(如圖3所示)：

圖3 Ransac算法流程

①根據(jù)式(2)計算最小循環(huán)次數(shù)N；

②計算擬合平面的標準差，如果其大于閾值δ0，重新確定平面的初值；否則繼續(xù)下一步；

③計算面內(nèi)點的數(shù)量，并與閾值點數(shù)比較，如大于閾值，則計算平面的標準差；否則返回上一步；

④重復(fù)步驟②、步驟③，直到根據(jù)判斷準則得出最佳平面；

⑤重復(fù)步驟①～步驟④，直到模型點的個數(shù)小于給定的閾值Nmin為止，得出最終的分割結(jié)果。

4 地質(zhì)產(chǎn)狀算法

本文采用空間解析幾何和向量代數(shù)法計算結(jié)構(gòu)面產(chǎn)狀，其原理為，過不在同一直線上三點確定的平面法向量的Z軸方位角為γ，其角度范圍在0°～360°取值，計算時不考慮γ＇的方向性，規(guī)定其取值范圍為0°～90°，由圖4可以證明：巖層的傾角δ可通過法向量與Z軸的夾角γ＇獲得(角度相等)，并且?guī)r層傾向可以由法向量在水平面(xoy)上的投影向量的方向算出[10]。

圖4 巖層產(chǎn)狀計算示意

4.1 面法向量算法原理

在獲取每個獨立結(jié)構(gòu)面點云數(shù)據(jù)后，本文采用最小二乘法對點云數(shù)據(jù)進行平面擬合[11]，本文假設(shè)誤差只存在于Z軸方向上，從而計算該平面的平面法向量，設(shè)點云擬合的平面方程為：

Z=aX+bY+c

(3)

當取Z為觀測值，n為觀測點數(shù)目時，可推算出其觀測值方程為：

Z+V=aX+bY+c

(4)

將上式改寫為誤差方程：

(5)

其中：

(6)

(7)

4.2 面產(chǎn)狀算法原理

(8)

規(guī)定δ的取值范圍為0°～90°，則：

(9)

(10)

由式(10)可簡化得：

(11)

根據(jù)下表對X和Y值的正負進行修正即可得巖層傾向θ，詳見表1[12]。

傾向基本角θ修正 表1

5 實例分析

本方法中涉及兩類算法，一是基于Ransac的點云結(jié)構(gòu)面分割算法，二是結(jié)構(gòu)面產(chǎn)狀算法，后者在計算過程中皆使用固定常量，利用計算機編程逐步計算即可。而前者Ransac算法中存在人為設(shè)定參數(shù)，直接影響到計算結(jié)果。故本章將依托龍慶峽景區(qū)邊坡防護工程獲取邊坡點云數(shù)據(jù)，以此算法進行測試驗證。龍慶峽距北京城區(qū) 85 km，主要為高山峽谷地貌，分布巖壁大都近直立，高約 200 m左右，節(jié)理發(fā)育明顯，適宜進行成果對比。

RANSAC的點云分割效果主要由兩閾值參數(shù)決定，分別為點到平面的距離閾值δ0和聚類面最小點云數(shù)量Nmin[13]，為取得最佳的分割結(jié)果，本文通過控制變量法試驗求取δ0與Nmin的最佳參數(shù)值。此試驗以龍慶峽邊坡峭壁點云數(shù)據(jù)為樣本，選擇 15 m×20 m的典型巖石出露面兩處，將距離閾值δ0設(shè)為 5 cm、7 cm、10 cm、12 cm、15 cm，將迭代條件中剩余點云數(shù)量Nmin設(shè)為全部點云的20%和固定數(shù)值100個兩種。比較分析分割效果得到統(tǒng)計圖如圖5、圖6所示。

圖5 點云分割運行耗時統(tǒng)計

圖6 點云分割節(jié)理面?zhèn)€數(shù)統(tǒng)計

軟件耗時方面，當距離閾值δ0越大，運行時間越短。精度方面，在δ0值不超過 12 cm時，大部分點云可以被識別，當δ0值達到 10 cm和Nmin設(shè)置為100個點云，點云分割效果達到最佳，如圖7所示。

圖7 點云分割結(jié)果

當結(jié)構(gòu)面分割完成，各結(jié)構(gòu)面點云文件可依次輸入產(chǎn)狀計算算法，并最終輸出其計算結(jié)果。為進一步評價此結(jié)果準確性，本文將此次結(jié)構(gòu)面產(chǎn)狀計算結(jié)果分組統(tǒng)計，并與人工現(xiàn)場地質(zhì)調(diào)查結(jié)果進行比對，對比結(jié)果如表2所示：

計算結(jié)果比對 表2

根據(jù)統(tǒng)計結(jié)果及現(xiàn)場調(diào)查成果，區(qū)內(nèi)發(fā)育的三組結(jié)構(gòu)面中第一組為巖層的層面，其他兩組為近垂直相交的陡傾節(jié)理，控制邊坡巖體結(jié)構(gòu)發(fā)育特征。統(tǒng)計結(jié)果與地質(zhì)調(diào)查數(shù)據(jù)基本一致。

6 結(jié) 語

本文提出的基于PCL的地質(zhì)信息提取方法，通過三維激光掃描儀所獲取的點云數(shù)據(jù)，經(jīng)過RANSAC算法、向量代數(shù)法，最終得到了巖體的相關(guān)地質(zhì)參數(shù)，此技術(shù)的出現(xiàn)為后續(xù)地質(zhì)調(diào)查與危巖體的治理和設(shè)計提供了可靠的資料。但本算法設(shè)置的閾值較多，而地質(zhì)巖體的產(chǎn)狀信息多與其特定的地理環(huán)境相關(guān)，故此算法需要在更多地區(qū)，更多巖性條件下進行測試，從而形成更加完善的閾值設(shè)定機制，并為下一步的三維地質(zhì)的分析與建模奠定堅實的基礎(chǔ)。