基于三維點(diǎn)云的結(jié)構(gòu)面產(chǎn)狀獲取方法研究

馮文凱，曾唯恐，程柯力，易小宇，焦隆新

(成都理工大學(xué) 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610059)

1 研究背景

巖體結(jié)構(gòu)面是巖體內(nèi)發(fā)育的具有不同形態(tài)的面狀地質(zhì)界面，對(duì)揭露地質(zhì)體結(jié)構(gòu)有重要意義。傳統(tǒng)上的調(diào)查方法需要人手持羅盤貼近巖體測(cè)量，以獲取其幾何要素。近年來，隨著近景攝影測(cè)量、三維激光掃描等技術(shù)的發(fā)展，實(shí)現(xiàn)了遠(yuǎn)程非接觸式地獲取結(jié)構(gòu)面點(diǎn)云數(shù)據(jù)。董秀軍等[1-4]利用三維激光掃描儀調(diào)查高陡邊坡，并基于PolyWorks軟件解譯結(jié)構(gòu)面，輔助完成地質(zhì)編錄和后期成圖。胡偉等[5]利用近景攝影測(cè)量技術(shù)試驗(yàn)前后結(jié)構(gòu)面三維形貌數(shù)據(jù)，分析影響結(jié)構(gòu)面剪切特性因素。周春霖等[6]基于雙目重構(gòu)技術(shù)對(duì)巖體結(jié)構(gòu)面進(jìn)行非接觸測(cè)量并得到結(jié)構(gòu)面產(chǎn)狀。王明常等[7]以VirtuoZo數(shù)字?jǐn)z影測(cè)量站為平臺(tái)提取結(jié)構(gòu)面特征點(diǎn)的三維坐標(biāo)，利用Matlab程序?qū)崿F(xiàn)巖體結(jié)構(gòu)面信息的快速獲取。目前對(duì)點(diǎn)云的操作往往要通過特殊的點(diǎn)云處理平臺(tái)，且這些平臺(tái)不具備針對(duì)地質(zhì)的專門解決方案，或者存在軟件版權(quán)的限制，阻礙了技術(shù)的推廣。

Python語(yǔ)言作為一種開源的程序設(shè)計(jì)語(yǔ)言，具備簡(jiǎn)潔性、可讀性、可拓展性等優(yōu)點(diǎn)。本文利用Python語(yǔ)言，并基于最小二乘法、主成分分析法的算法原理[8-10]，構(gòu)建了一套擬合平面的方法，實(shí)現(xiàn)了點(diǎn)云平面坐標(biāo)向結(jié)構(gòu)面產(chǎn)狀的計(jì)算，為結(jié)構(gòu)面點(diǎn)云的解譯提供了新的輕量化方案。

2 基于Python的結(jié)構(gòu)面產(chǎn)狀獲取流程

利用近景攝影測(cè)量、三維激光掃描獲取到的結(jié)構(gòu)面數(shù)據(jù)通常稱為點(diǎn)云，它是一系列離散三維點(diǎn)的集合，共同組成了結(jié)構(gòu)面的三維輪廓。為了通過點(diǎn)云數(shù)據(jù)獲取產(chǎn)狀信息，首先需要擬合出一個(gè)虛擬平面以得到平面方程。

2.1 最小二乘法擬合平面

2.1.1 最小二乘法擬合平面原理

最小二乘法最早出自法國(guó)數(shù)學(xué)家阿德里安，他提出讓誤差的平方和最小的y值就是真值(圖1)。因此，擬合平面的過程可以視作尋求最小化誤差下的平面函數(shù)。對(duì)于組成結(jié)構(gòu)面的起伏的點(diǎn)云，每一個(gè)點(diǎn)并不總是精準(zhǔn)地落在擬合平面中，我們用點(diǎn)到擬合平面的距離在Z軸方向上的投影來衡量誤差，那么滿足該值最小的平面即為所需要的擬合平面。

圖1 最小二乘法原理Fig.1 Principle of ordinary least square

2.1.2 求解思路

設(shè)點(diǎn)云數(shù)據(jù)集Q=(xi，yi，zi)，(i=1，2，3，…，n)，其中x，y，z為3個(gè)維度上的位置坐標(biāo)。設(shè)要求的平面方程為

z=ax+by+c。

(1)

式中a、b、c均為平面方程參數(shù)。

根據(jù)最小二乘法定律，各點(diǎn)到平面的距離平方和S應(yīng)為最小，即

(2)

通過求其偏導(dǎo)并設(shè)為0以求參數(shù)，得如下方程組：

(3)

移項(xiàng)后用矩陣形式表達(dá)為

式(4)可轉(zhuǎn)換為

AX=b。

(5)

在式子兩邊左乘A的逆矩陣A-1，則有

X=A-1b。

(6)

從而解出平面參數(shù)X(a，b，c)。

2.1.3 Python程序流程

本文選擇Python語(yǔ)言設(shè)計(jì)算法，Python具備簡(jiǎn)潔性、可讀性、可拓展性等優(yōu)點(diǎn)，通過Numpy科學(xué)計(jì)算拓展庫(kù)，使其具備數(shù)值計(jì)算的功能。計(jì)算步驟如下：

(1)導(dǎo)入點(diǎn)云數(shù)據(jù)集Q=(xi，yi，zi)，(i=1，2，3，…，n)。

(2)創(chuàng)建其系數(shù)矩陣A、常數(shù)列b，其表達(dá)式為：

(7)

(8)

(3)利用Python中np.linage.linv()求逆函數(shù)計(jì)算A的逆矩陣，即

A.inv=np.linalg.inv(A) 。

(9)

(4)用A.inv點(diǎn)乘b，解得平面參數(shù)X，其表達(dá)式為

X=np.dot(A.inv,b)=(a,b,c) 。

(10)

2.2 主成分分析法擬合平面

2.2.1 主成分分析法擬合平面原理

主成分分析法(Principal Component Analysis，PCA)的核心思想是將一組變量線型變換為另一組不相關(guān)的正交變量，使得變換后的第一變量方差最大，其余變量方差依次遞減至可以忽略。由于方差代表了信息量，也就意味著第一變量所包含的信息量最大，前幾組變量已經(jīng)包含了絕大部分信息量，足以替換原始變量，最終達(dá)到降維的目的。

從幾何角度來說，PCA對(duì)點(diǎn)云數(shù)據(jù)的降維可看作是將點(diǎn)云投影到另一個(gè)三維直角坐標(biāo)系(圖2)，3個(gè)新的坐標(biāo)軸代表了3個(gè)最大主成分。如圖2所示，在新的坐標(biāo)系下，點(diǎn)云在X′O′Y′平面上方差值最大，在Z軸上方差最小，X′O′Y′平面即為最佳擬合平面，Z軸方向即為平面法向量。

圖2 PCA變換原理Fig.2 Principle of princi- pal components analysis

2.2.2 求解思路

對(duì)n個(gè)包含三維坐標(biāo)信息的點(diǎn)云數(shù)據(jù)可以用n×3的矩陣X來表示，即

(11)

為了求其主成分方向，需要找到樣本協(xié)方差矩陣XTX的最大3個(gè)特征向量。由于特征向量具有兩兩相互正交的特征，因此前2個(gè)最大的特征向量代表了最大主成分，二者組成了擬合平面。

Python提供有“np.linalg.eig()”的函數(shù)求解特征向量，但對(duì)高維數(shù)據(jù)的暴力求解會(huì)導(dǎo)致計(jì)算量偏大。而線性代數(shù)中的奇異值分解(Singular Value Decomposition，SVD)法也能求解特征向量，且不需要先計(jì)算協(xié)方差矩陣XTX，因此我們引入SVD法。

2.2.3 python程序流程

(1)導(dǎo)入n×3階點(diǎn)云矩陣Xn×3并做中心化處理。

(2)利用np.linage.svd()函數(shù)對(duì)X作奇異值分解，即

X.svd=np.linalg.svd(X) 。

(12)

得到的X.svd對(duì)象包含U，Σ，VT三個(gè)矩陣，VT稱為右奇異矩陣，由樣本的特征向量所組成,即

X.svd=UΣVT。

(13)

(3)取右奇異矩陣VT為

VT=X.svd[2] 。

(14)

(4)取右奇異矩陣的第3行向量(第三主成分)：

n=VT[2] 。

(15)

根據(jù)前文介紹，點(diǎn)云矩陣的第三主成分即為平面的法向量n，通過平面法向量可反求得平面方程。

2.3 試驗(yàn)點(diǎn)檢驗(yàn)算法

在第2.2節(jié)中，利用Python語(yǔ)言分別構(gòu)建了最小二乘法、主成分分析法(SVD分解)2種平面擬合算法。為了對(duì)比2種算法試驗(yàn)效果，生成了1 000個(gè)虛擬三維點(diǎn)作為試驗(yàn)集，試驗(yàn)點(diǎn)均滿足平面方程：2x+3y+5=z，其平面參數(shù)為(2，3，5)。為了模擬出結(jié)構(gòu)面本身的凹凸不平，同時(shí)加入了100個(gè)與平面稍有偏移并呈正態(tài)分布的噪聲點(diǎn)(圖3)。

圖3 建立試驗(yàn)點(diǎn)點(diǎn)云Fig.3 Establishment of test point cloud dataset

用上述試驗(yàn)點(diǎn)擬合平面，擬合結(jié)果如圖4所示，采用最小二乘法得到的平面參數(shù)為(2.047 136，3.053 119，4.962 118)，采用主成分分析法得到的平面參數(shù)為(2.007 714，2.991 216，4.995 733)。

圖4 兩種算法擬合結(jié)果Fig.4 Fitting results of two algorithms

由圖4可知，利用主成分分析法得到的平面方程更加精確，原因在于在最小二乘法中，需要我們不斷優(yōu)化的誤差值d1，計(jì)算的是點(diǎn)與擬合平面在Z軸方向上的投影距離(圖5)，而實(shí)際該值用沿平面法向量上的距離差值d2來衡量更加準(zhǔn)確。主成分分析法則規(guī)避了這一點(diǎn)。

圖5 誤差來源Fig.5 Source of error

2.4 轉(zhuǎn)化為結(jié)構(gòu)面產(chǎn)狀

為了滿足地質(zhì)工作要求，需要將前文所獲得的平面方程參數(shù)轉(zhuǎn)化為結(jié)構(gòu)面產(chǎn)狀。

一般來說，X軸代表了正東方向，Y軸代表了正北方向，Z軸指向垂直向上，它們?cè)诳臻g中的對(duì)應(yīng)關(guān)系見圖6，圖中n(nx，ny，nz)為平面的單位法向量。

圖6 空間大地直角坐標(biāo)系Fig.6 Geodetic coordinate system

根據(jù)結(jié)構(gòu)面產(chǎn)狀的定義和圖中所示幾何關(guān)系可知，傾角α是結(jié)構(gòu)面與水平面間夾角，也即單位法向量n在Z軸上的分量nz與Z軸的夾角。計(jì)算公式為

α=arccos(nz) 。

(16)

傾向β是其單位法向量n在水平方向上的投影與正北方向之間的夾角，根據(jù)n投影在XOY平面上象限的不同有4種情況，其計(jì)算公式可化簡(jiǎn)為如下2種：

當(dāng)nx≥0時(shí)，

(17)

當(dāng)nx≤0時(shí)，

(18)

由上述計(jì)算過程可知，傾向傾角的求解均需要單位法向量n，n按式(19)求得，其中a，b，c為平面方程(式1)的平面參數(shù)。

n(nx,ny,nz)=

(19)

通過上述理論在Python中編程實(shí)現(xiàn)了從平面方程向產(chǎn)狀的轉(zhuǎn)換，從而完成了結(jié)構(gòu)面點(diǎn)云擬合平面并計(jì)算產(chǎn)狀的完整過程。

3 實(shí)例驗(yàn)證

通過國(guó)際公開實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證這套算法的適用性，結(jié)構(gòu)面點(diǎn)云數(shù)據(jù)采集于美國(guó)科羅拉多州烏雷地區(qū)一處裸露基巖[11](圖7)。

圖7 結(jié)構(gòu)面現(xiàn)場(chǎng)照片F(xiàn)ig.7 Photo of structural planes

3.1 截取點(diǎn)云代入計(jì)算

如圖7所示，J1—J5分別代表了5塊產(chǎn)狀各異的結(jié)構(gòu)面。運(yùn)用“cloudcompare”軟件平臺(tái)從點(diǎn)云數(shù)據(jù)中分別截取對(duì)應(yīng)的結(jié)構(gòu)面點(diǎn)云，導(dǎo)入前文程序中進(jìn)行計(jì)算，計(jì)算結(jié)果與參考文獻(xiàn)[11]中結(jié)果的對(duì)比見表1。

表1 產(chǎn)狀結(jié)果對(duì)比Table 1 Comparison of occurrence results

從表1可以看出：對(duì)于規(guī)則的結(jié)構(gòu)面點(diǎn)云數(shù)據(jù)，計(jì)算結(jié)果與實(shí)際數(shù)據(jù)差距<2°，2種方法都有較好的擬合效果。同時(shí)觀察到，二者之間的差別基本上<1°，但對(duì)于傾向較小的緩傾層面，該差值能縮小至0.1°，并隨著傾向的變大而逐漸增加。其原因如第2.3節(jié)所述，是衡量誤差的距離值取法不同所致，該問題在平緩結(jié)構(gòu)面上表現(xiàn)得不明顯，但差異尚不足以影響工程實(shí)際應(yīng)用。

3.2 結(jié)合聚類算法自動(dòng)計(jì)算

上述算法的實(shí)施需要事先設(shè)置好結(jié)構(gòu)面點(diǎn)云數(shù)據(jù)，工程人員可以根據(jù)自己主觀判斷計(jì)算任意一塊結(jié)構(gòu)面，也可以只提取結(jié)構(gòu)面中較為平整的一小塊，使得結(jié)果更加精準(zhǔn)。若是采用無監(jiān)督聚類算法提前對(duì)結(jié)構(gòu)面自動(dòng)聚類既可以排除該過程的人工干預(yù)，同時(shí)也能提高整體解譯效率。

本文選取譜聚類算法對(duì)點(diǎn)云進(jìn)行無監(jiān)督分類。該算法同樣基于Python實(shí)現(xiàn)，這是一種基于圖論的聚類方法[12]，能識(shí)別任意形狀的樣本空間且能獲得全局最優(yōu)解，相比于K-Means等常用的聚類算法，能有效應(yīng)用于三維點(diǎn)云此類非凸分布樣本。計(jì)算過程需要首先去除點(diǎn)云數(shù)據(jù)中無關(guān)部分，再利用譜聚類算法分簇，結(jié)果如圖8(坐標(biāo)軸為3個(gè)正交方向)所示，除了J1所在平面不連續(xù)而分成了4塊后，其余結(jié)構(gòu)面都被較好地識(shí)別，呈現(xiàn)較好的分割效果。

圖8 譜聚類分簇結(jié)果Fig.8 Clustering results by spectral clustering

對(duì)各組結(jié)構(gòu)面進(jìn)行循環(huán)計(jì)算后，得到了2種算法的計(jì)算結(jié)果，數(shù)據(jù)對(duì)比見表2。

觀察表2可以看出，無監(jiān)督聚類后計(jì)算的結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)較真實(shí)數(shù)據(jù)出現(xiàn)一定偏差，但誤差基本<5°，其中J5結(jié)構(gòu)面本身起伏較大，難以估計(jì)，故偏差較大。在此方法中，聚類好壞直接影響到結(jié)果的準(zhǔn)確性，同時(shí)也能看出，在有干擾點(diǎn)被錯(cuò)誤聚類到結(jié)構(gòu)面數(shù)據(jù)中時(shí)，主成分分析法較最小二乘法更加穩(wěn)定。

表2 聚類分簇后產(chǎn)狀結(jié)果對(duì)比Table 2 Comparison of occurrence results with clustering results

4 結(jié) 論

本文利用Python語(yǔ)言分別構(gòu)建了傳統(tǒng)最小二乘法和主成分分析法(SVD分解)2種算法以計(jì)算點(diǎn)云結(jié)構(gòu)面產(chǎn)狀。在實(shí)例中進(jìn)行運(yùn)用并得到以下結(jié)論：

(1)2種方法對(duì)于出露條件較好、起伏較小的結(jié)構(gòu)面點(diǎn)云數(shù)據(jù)均有較好的擬合效果，在點(diǎn)云輔助地質(zhì)調(diào)查中具有廣泛應(yīng)用前景。受誤差來源影響，主成分分析法要整體優(yōu)于最小二乘法。

(2)由于上述擬合方法的假設(shè)前提是待檢測(cè)的結(jié)構(gòu)面點(diǎn)云本身就是標(biāo)準(zhǔn)的平面，因此當(dāng)點(diǎn)云數(shù)據(jù)越趨近于完整平面時(shí)，計(jì)算結(jié)果越好。對(duì)于出露條件不好或是掃描數(shù)據(jù)失真的點(diǎn)云數(shù)據(jù)，通過該方法計(jì)算得到的結(jié)果則偏差較大，應(yīng)及時(shí)人工輔助判斷。

(3)應(yīng)用無監(jiān)督聚類算法自動(dòng)給結(jié)構(gòu)面分組能有效提高產(chǎn)狀解譯速度，減少人為干預(yù)工作量，但會(huì)增大計(jì)算結(jié)果的誤差。提高聚類算法的準(zhǔn)確度能極大優(yōu)化數(shù)據(jù)結(jié)果，這同時(shí)也是未來針對(duì)點(diǎn)云結(jié)構(gòu)面的重要研究方向。