三維巖體點云數據地質信息提取方法

王喆，劉鵬飛，王凱

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

1 引 言

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

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

利用三維激光掃描技術可實現無接觸地直接采集目標巖體的點云數據，數據經過處理后得到采集目標的三維點云模型。由于該模型真實地逼近于現場情況，技術人員可在后期處理軟件中對危巖體幾何參數進行量測和計算，得到巖體結構面產狀，從而完成一系列地質分析。故此方法在地質勘察工作中，尤其是大高差山區危巖地區，可基本代替傳統人工勘測方法，對地質體信息提取工作具有極大的價值[3]。

2 PCL簡介

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

圖1 PCL架構圖

2.1 1PCD文件

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

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

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

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

2.2 PCD文件操作

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

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

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

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

圖2 類File Reader和File-Writer的繼承關系

3 平面分割算法

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

P=1-(1-εk)n

(1)

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

(2)

算法迭代的核心邏輯為：

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

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

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

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

圖3 Ransac算法流程

①根據式(2)計算最小循環次數N；

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

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

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

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

4 地質產狀算法

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

圖4 巖層產狀計算示意

4.1 面法向量算法原理

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

Z=aX+bY+c

(3)

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

Z+V=aX+bY+c

(4)

將上式改寫為誤差方程：

(5)

其中：

(6)

(7)

4.2 面產狀算法原理

(8)

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

(9)

(10)

由式(10)可簡化得：

(11)

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

傾向基本角θ修正 表1

5 實例分析

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

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

圖5 點云分割運行耗時統計

圖6 點云分割節理面個數統計

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

圖7 點云分割結果

當結構面分割完成，各結構面點云文件可依次輸入產狀計算算法，并最終輸出其計算結果。為進一步評價此結果準確性，本文將此次結構面產狀計算結果分組統計，并與人工現場地質調查結果進行比對，對比結果如表2所示：

計算結果比對 表2

根據統計結果及現場調查成果，區內發育的三組結構面中第一組為巖層的層面，其他兩組為近垂直相交的陡傾節理，控制邊坡巖體結構發育特征。統計結果與地質調查數據基本一致。

6 結 語

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