三維巖體點(diǎn)云數(shù)據(jù)地質(zhì)信息提取方法

王喆，劉鵬飛，王凱

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

1 引 言

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

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

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

2 PCL簡(jiǎn)介

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

圖1 PCL架構(gòu)圖

2.1 1PCD文件

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

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

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

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

2.2 PCD文件操作

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

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

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

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

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

3 平面分割算法

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

P=1-(1-εk)n

(1)

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

(2)

算法迭代的核心邏輯為：

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

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

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

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

圖3 Ransac算法流程

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

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

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

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

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

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

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

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

4.1 面法向量算法原理

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

Z=aX+bY+c

(3)

當(dāng)取Z為觀測(cè)值，n為觀測(cè)點(diǎn)數(shù)目時(shí)，可推算出其觀測(cè)值方程為：

Z+V=aX+bY+c

(4)

將上式改寫(xiě)為誤差方程：

(5)

其中：

(6)

(7)

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

(8)

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

(9)

(10)

由式(10)可簡(jiǎn)化得：

(11)

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

傾向基本角θ修正 表1

5 實(shí)例分析

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

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

圖5 點(diǎn)云分割運(yùn)行耗時(shí)統(tǒng)計(jì)

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

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

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

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

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

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

6 結(jié) 語(yǔ)

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