高陡邊坡巖體結(jié)構(gòu)信息無人機(jī)識(shí)別與離散元數(shù)值模擬*

沙 鵬 孔德珩 王紹亮 伍法權(quán)

(①紹興文理學(xué)院，浙江省巖石力學(xué)與地質(zhì)災(zāi)害重點(diǎn)實(shí)驗(yàn)室，紹興 312000，中國)(②同濟(jì)大學(xué)，土木工程學(xué)院，上海 200092，中國)(③浙江省工程勘察設(shè)計(jì)院集團(tuán)有限公司，寧波 315012，中國)

0 引 言

巖質(zhì)高陡邊坡的穩(wěn)定性評價(jià)需要地質(zhì)工作者在現(xiàn)場調(diào)查、測量、分析的基礎(chǔ)上，明確邊坡巖體結(jié)構(gòu)的空間展布狀態(tài)(孫廣忠，1993；黃潤秋等，2004)。對于復(fù)雜艱險(xiǎn)山區(qū)的自然高邊坡，測量人員根本無法到達(dá)所需要測量的區(qū)域，傳統(tǒng)測量方法無法滿足大斷面巖體結(jié)構(gòu)測量(于天亮，2009)。對于水電、礦山等工程開挖形成的高陡邊坡，現(xiàn)場要求快速完成地質(zhì)編錄工作，不僅需要投入大量人力，而且崩塌等偶發(fā)情況可能對測量人員帶來安全隱患。因此，巖體工程領(lǐng)域需要引進(jìn)快速、高精度、高效率且對地形地貌有較好適應(yīng)性的調(diào)查技術(shù)。

隨著三維空間影像技術(shù)的不斷成熟與完善，如數(shù)字近景攝影測量(王鳳艷等，2008)、無人機(jī)低空遙感技術(shù)(Giordan et al.，2018；張愷等，2019)和三維激光掃描技術(shù)(Jaboyedoff et al.，2012；Riquelme et al.，2015；梁玉飛等，2021)等，為解決上述問題提供了新的思路。其中：無人機(jī)低空遙感技術(shù)和三維激光掃描技術(shù)能夠遠(yuǎn)距離、高精度快速獲取點(diǎn)云數(shù)據(jù)，便于進(jìn)行三維重建和算法智能解譯，這為實(shí)現(xiàn)結(jié)構(gòu)面的精細(xì)化識(shí)別和參數(shù)提取提供了技術(shù)支持(Feng et al.，2014；Chen et al.，2016)。基于三維點(diǎn)云模型識(shí)別巖體結(jié)構(gòu)面的研究日益增多。基于不規(guī)則三角網(wǎng)模型以及不同的幾何特征聚類分析，如模糊K-均值算法與改進(jìn)K-均值算法等實(shí)現(xiàn)結(jié)構(gòu)面產(chǎn)狀與分組的自動(dòng)識(shí)別(Slob，2010；V?ge et al.，2013；Liu L et al.，2019)。Kong et al.(2020)利用通過快速搜索和查找密度峰(CFSFDP)算法進(jìn)行聚類來自動(dòng)檢測聚類的數(shù)量，提高了準(zhǔn)確性并縮短計(jì)算時(shí)間。結(jié)構(gòu)面的間距和線密度可通過虛擬精測線法、平面間距離計(jì)算等方法在三維點(diǎn)云數(shù)據(jù)中提取(葛云峰等，2017)。可以看出，大部分前人研究主要以某個(gè)結(jié)構(gòu)參數(shù)，如產(chǎn)狀、間距等參數(shù)的識(shí)別方法為主，尚未完整建立巖體結(jié)構(gòu)模型。

另一方面，如何建立精確的數(shù)值計(jì)算模型開展數(shù)值模擬也是邊坡穩(wěn)定性分析中的研究熱點(diǎn)。不少學(xué)者通過致力于進(jìn)行三維地學(xué)建模與數(shù)值模擬的耦合，彌補(bǔ)數(shù)值模型在建立復(fù)雜三維地質(zhì)模型上的不便(崔芳鵬等，2008；徐文杰等，2008)。三維點(diǎn)云數(shù)據(jù)可以獲取邊坡地形與關(guān)鍵塊體信息，通過相關(guān)軟件將DEM導(dǎo)入數(shù)值模擬軟件，分析不同工況下的邊坡形變或不穩(wěn)定塊體運(yùn)動(dòng)，實(shí)現(xiàn)邊坡穩(wěn)定性評價(jià)的目的(Havaej et al.，2016；Liu C et al.，2019；金愛兵等，2021)，但上述研究成果尚未考慮結(jié)構(gòu)面對巖質(zhì)邊坡穩(wěn)定性的影響。如何利用三維點(diǎn)云數(shù)據(jù)建立巖體結(jié)構(gòu)模型，結(jié)合DEM快速建立三維離散元數(shù)值分析模型需要進(jìn)一步開展研究。基于此，本文結(jié)合露天礦山高邊坡無人機(jī)地質(zhì)調(diào)查，基于獲取的三維點(diǎn)云數(shù)據(jù)實(shí)現(xiàn)巖體結(jié)構(gòu)幾何信息的自動(dòng)識(shí)別。采用曲面重構(gòu)和“Rhino-Griddle-3DEC”聯(lián)合建模方法，建立邊坡三維離散元數(shù)值模型。結(jié)合巖體結(jié)構(gòu)統(tǒng)計(jì)模型，構(gòu)建離散裂隙網(wǎng)絡(luò)(DFN)切割邊坡模型，開展邊坡穩(wěn)定性分析計(jì)算。

1 工程背景

圖1 銀山銅礦露天邊坡全景與巖體分區(qū)Fig.1 Panorama of Yinshan copper open pit slope and distribution of rock massesa.高陡巖質(zhì)邊坡；b.各分區(qū)巖體分布

邊坡巖體主要為千枚巖，其次有安斑巖、石英斑巖、石英閃長巖等，具體巖體分布特征如圖1b。其中：I區(qū)-Ⅱ區(qū)的巖體主要為絹云母千枚巖，片理發(fā)育程度不高，構(gòu)造裂隙不甚發(fā)育。Ⅲ區(qū)千枚巖片理發(fā)育，巖體較破碎，構(gòu)造裂隙較發(fā)育，且根據(jù)英安斑巖、石英斑巖等分布情況繼續(xù)劃分成A、B、C 3部分。Ⅳ區(qū)內(nèi)千枚巖發(fā)育一組片理，主要產(chǎn)狀為285°∠82°。

表 1 大疆Phantom 4 Pro 四旋翼無人機(jī)主要參數(shù)Table1 Main parameters of DJI Phantom 4 Pro UAV system

2 無人機(jī)點(diǎn)云數(shù)據(jù)處理

2.1 測試方案

圖2 礦山邊坡無人機(jī)測量Fig.2 UAV measurements of the open pit slope

2.2 無人機(jī)點(diǎn)云數(shù)據(jù)獲取及處理

采用運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)法(Structure-from-Motion，SfM)，篩選無人機(jī)航拍圖像導(dǎo)入Context Capture軟件，生成礦山邊坡各個(gè)分區(qū)的高精度三維點(diǎn)云模型，如圖3所示。SfM法是通過圖像集中的同名點(diǎn)來估計(jì)靜止場景中運(yùn)動(dòng)相機(jī)的相對參數(shù)，使用相機(jī)參數(shù)以及同名點(diǎn)間對極幾何關(guān)系來獲取三維點(diǎn)云模型，恢復(fù)三維場景的一種方法(Westoby et al.，2012)。這個(gè)過程涉及三維幾何和攝像機(jī)運(yùn)動(dòng)姿態(tài)的同時(shí)估計(jì)運(yùn)算，因此稱為運(yùn)動(dòng)恢復(fù)結(jié)構(gòu)(孫鈺杰，2016)。

圖3 露天礦山邊坡各區(qū)無人機(jī)點(diǎn)云圖像Fig.3 UAV photogrammetry of districts of open-pit slope

由于分區(qū)較多，受篇幅限制下文只針對北東邊坡Ⅲ-B區(qū)域進(jìn)行詳細(xì)分析。在Ⅲ-B區(qū)域(東處高邊坡120m平臺(tái))共設(shè)地面控制標(biāo)靶站點(diǎn)4處，獲得航片872張。為了獲得精確的位置，測量數(shù)據(jù)需要利用機(jī)載鏡頭的投影矩陣來計(jì)算出匹配點(diǎn)的三維坐標(biāo)。此方法流程圖如圖4所示，包括數(shù)據(jù)預(yù)處理、特征點(diǎn)提取匹配、稀疏重構(gòu)、密集重構(gòu)。

圖4 SfM技術(shù)流程圖Fig.4 The flow chart of SfM method

(1)數(shù)據(jù)預(yù)處理：完成相機(jī)的內(nèi)部參數(shù)標(biāo)定和檢驗(yàn)拍攝照片的質(zhì)量。

(2)特征點(diǎn)提取匹配：特征點(diǎn)匹配是尋找不同影像中相同特征點(diǎn)的過程。考慮到使用無人機(jī)進(jìn)行拍攝，會(huì)出現(xiàn)照片尺度、旋轉(zhuǎn)角度、環(huán)境光強(qiáng)度、模糊程度的變化，因此采用尺度不變特征變換算法(SIFT，Lowe，2004)對航攝圖像進(jìn)行特征提取匹配。該算法提取特征點(diǎn)時(shí)具有圖像旋轉(zhuǎn)、平移、縮放以及仿射不變型和一定的抗光照強(qiáng)度及視角變化能力，具體操作步驟為：

1)使用高斯卷積函數(shù)對無人機(jī)航拍圖像進(jìn)行尺度變換，獲得圖像在不同尺度下的表達(dá)序列。

2)利用尺度空間高斯差分的方程形式對圖像進(jìn)行卷積求取極值，這些極值所對應(yīng)的極值點(diǎn)即是圖像的特征點(diǎn)。

3)將坐標(biāo)軸旋轉(zhuǎn)到特征點(diǎn)方向，以保證旋轉(zhuǎn)角度不變，然后計(jì)算SIFT特征向量。

4)結(jié)合無人機(jī)采集圖像時(shí)的RTK-GPS坐標(biāo)信息、慣性測量單元位姿角度信息建立圖像間的拓?fù)浣Y(jié)構(gòu)，利用最鄰近算法計(jì)算出圖像間特征點(diǎn)的對應(yīng)關(guān)系完成匹配工作。

(3)稀疏重構(gòu)：利用上述匹配出具有同名特征的點(diǎn)集對測量場景進(jìn)行增量式重構(gòu)，得到初步的三維點(diǎn)云模型。基本流程是：

1)首先求取內(nèi)外參數(shù)，整個(gè)測量過程中，無人機(jī)通過同一鏡頭拍攝，因此圖像的內(nèi)部參數(shù)是固定的，只需計(jì)算獲得圖像的外部參數(shù)即圖像間的旋轉(zhuǎn)矩陣和平移向量。

2)根據(jù)相鄰圖片的內(nèi)外參數(shù)，采用三角定位方法解算空間點(diǎn)的三維坐標(biāo)。

3)由于特征點(diǎn)匹配并不是絕對精確，因此步驟②中求得三維點(diǎn)坐標(biāo)存在誤差。采用光束平差法，通過逐步迭代最小化投影點(diǎn)和觀測圖像點(diǎn)之間的重投影誤差，可以得到機(jī)載鏡頭最佳姿態(tài)和所測場景的三維點(diǎn)云坐標(biāo)。

4)密集重構(gòu)：由于稀疏重構(gòu)所恢復(fù)的三維點(diǎn)云模型無法對巖體這類不規(guī)則物體進(jìn)行精細(xì)描述。因此，需要對點(diǎn)云模型進(jìn)行密集重構(gòu)，即增加點(diǎn)云密度提高模型的精細(xì)程度。采用Furukawa et al.(2010)提出的稠密重建算法(CMVS-PMVS)，采用聚簇分類去除原始圖像組中的冗余圖像；選取可信度較高的稀疏匹配點(diǎn)生成一系列稀疏面片；將對應(yīng)的圖像區(qū)域進(jìn)行反復(fù)的擴(kuò)散、濾波，最終得到更為精細(xì)的密集點(diǎn)云數(shù)據(jù)。

如圖5所示，通過上述4步處理，Ⅲ-B區(qū)120m平臺(tái)共生成點(diǎn)云數(shù)量64，353，602個(gè)，構(gòu)建了長度約20m，高度約12m的三維點(diǎn)云模型。截取長18m、高6m且結(jié)構(gòu)面出露較好的矩形區(qū)域，利用點(diǎn)云數(shù)據(jù)包含的真實(shí)空間坐標(biāo)信息開展后續(xù)的結(jié)構(gòu)面的智能識(shí)別和幾何信息解譯?；舅悸分饕峭ㄟ^算法自動(dòng)處理整體點(diǎn)云數(shù)據(jù)，篩選出代表結(jié)構(gòu)面的點(diǎn)，實(shí)現(xiàn)程序化識(shí)別結(jié)構(gòu)面。

圖5 Ⅲ-B區(qū)部分航拍照片(a)；Ⅲ-B區(qū)局部點(diǎn)云(b)Fig.5 Partial aerial photos of area Ⅲ-B(a)；Partial point cloud of Ⅲ-B area(b)

3 結(jié)構(gòu)面自動(dòng)提取方法

3.1 點(diǎn)法向量計(jì)算

本文采用基于Riquelme et al.(2014)提出的迭代加權(quán)平面擬合(IRPF)方法，計(jì)算區(qū)域露頭點(diǎn)云集合內(nèi)的各點(diǎn)法向量。具體可分為以下3個(gè)部分：

(1)估算每個(gè)點(diǎn)的初始法向量。利用式(1)進(jìn)行迭代的加權(quán)平面擬合，其中將σd設(shè)置為局部點(diǎn)間距δl，將σr設(shè)置為最大擬合殘差的1/3，每個(gè)點(diǎn)的初始擬合剩余權(quán)重設(shè)置為1。當(dāng)涉及的權(quán)重穩(wěn)定或迭代次數(shù)超過所設(shè)限制的20時(shí)，迭代即終止。

(1)

(2)

(2)估算局部自適應(yīng)帶寬。對于任意點(diǎn)，將其鄰近點(diǎn)的法向量在3δl之間的距離內(nèi)進(jìn)行聚類，以便自動(dòng)確定其特征系數(shù)。生成尖銳的特征點(diǎn)百分比曲線時(shí)，角度閾值從0開始，以2.5°為增量增加。找到合適的閾值并估計(jì)每個(gè)點(diǎn)的特征系數(shù)后，根據(jù)式(3)和式(4)設(shè)置3個(gè)權(quán)重帶寬。

(3)

(4)

因此，對于平滑區(qū)域中的點(diǎn)(Fcoe=0)，距離帶寬為δl。對于尖銳特征點(diǎn)(Fcoe=1)，距離帶寬2δl將包含更多鄰近點(diǎn)，以補(bǔ)償由于剩余和正常差值權(quán)重而丟棄的鄰近點(diǎn)的缺失。

(3)改善初始法向量。采用式(5)來改善初始法向量，迭代計(jì)算通過與式(1)中相同的方式終止。

(5)

3.2 識(shí)別同組結(jié)構(gòu)面并聚類分割

根據(jù)求得的各點(diǎn)法向量對點(diǎn)云數(shù)據(jù)進(jìn)行相似性判斷和聚類，確定主要點(diǎn)云集合即為同組結(jié)構(gòu)面，并在赤平投影圖中進(jìn)行統(tǒng)計(jì)分析，由此得到點(diǎn)云集合的主要方向即該組結(jié)構(gòu)面的優(yōu)勢產(chǎn)狀，整個(gè)過程采用基于劃分的K-Means算法(Macqueen，1967)。通過對點(diǎn)云集合進(jìn)行聚類分割與平面擬合，得到的最佳擬合平面即為結(jié)構(gòu)面，并可通過法向量計(jì)算對應(yīng)的結(jié)構(gòu)面產(chǎn)狀。

為了識(shí)別單條結(jié)構(gòu)面，需要從點(diǎn)云集合中識(shí)別出結(jié)構(gòu)面所在平面的空間點(diǎn)簇。本文采用基于密度梯度的聚類算法(Wang et al.，2013)，從上述聚類后的點(diǎn)云集合進(jìn)行分割提取。具體步驟如下：

(1)初始化，計(jì)算密度分布：計(jì)算各點(diǎn)的d近鄰，以其d近鄰的平均距離作為該點(diǎn)的密度。

(2)獲取不動(dòng)點(diǎn)，獲得原始聚類：隨機(jī)選擇某一未分類的點(diǎn)O，獲取其d近鄰中密度最大的點(diǎn)P，根據(jù)以下兩種情況比較點(diǎn)P與點(diǎn)O的密度：a.當(dāng)點(diǎn)P密度小于點(diǎn)O密度，即D(P，d)的數(shù)值大于D(O，d)的數(shù)值時(shí)，則點(diǎn)O為一不動(dòng)點(diǎn)，并賦予新的類別號(hào)。b.當(dāng)點(diǎn)P密度大于等于點(diǎn)O密度時(shí)，根據(jù)點(diǎn)P的聚類情況進(jìn)行判斷，若該點(diǎn)類別號(hào)已經(jīng)確定，可將此聚類號(hào)作為該點(diǎn)的聚類號(hào)；若該點(diǎn)類別號(hào)沒有確定，點(diǎn)P可看作起始點(diǎn)的不動(dòng)點(diǎn)。確定相應(yīng)的不動(dòng)點(diǎn)后，將該不動(dòng)點(diǎn)的類別號(hào)賦予該路徑下的所有點(diǎn)，由此繼續(xù)循環(huán)至所有點(diǎn)都進(jìn)行分類。

(3)調(diào)整邊界點(diǎn)：選擇邊界點(diǎn)的d個(gè)近鄰點(diǎn)的類別號(hào)作為該邊界點(diǎn)的類別號(hào)。

(4)計(jì)算邊界點(diǎn)對各聚類的分布情況進(jìn)行合并及分割：將單條結(jié)構(gòu)面所在平面的空間點(diǎn)簇分割提取出來，接下去只需將這些點(diǎn)簇進(jìn)行平面擬合，即為所需結(jié)構(gòu)面。采用最小二乘法擬合平面獲得結(jié)構(gòu)面所在平面的方程：

Ax+By+Cz+D=0

(6)

(7)

如果記：

(8)

那么：

z=a0x+a1y+a2

(9)

式中：a0，a1，a2為方程系數(shù)。對于屬于一系列的n個(gè)點(diǎn)的坐標(biāo)(xi，yi，zi)，i=1，2，3…n，使用最小二乘法擬合平面，則滿足各點(diǎn)到擬合平面的距離平方和S為最小值，S的計(jì)算公式為：

(10)

若要使S為最小值，應(yīng)滿足下式：

(11)

聯(lián)立可得：

(12)

解得：

(13)

通過代入a0，a1，a2即可得到單位法向量的分量。

對結(jié)構(gòu)面進(jìn)行平面擬合后，即可通過計(jì)算結(jié)構(gòu)面的單位的法向量以得到結(jié)構(gòu)面的產(chǎn)狀：

(14)

其中：ua，ub，uc為擬合平面單位法向量的分量。

3.3 幾何參數(shù)獲取

幾何參數(shù)主要是結(jié)構(gòu)面位置、跡長和間距。其中結(jié)構(gòu)面的位置由同組點(diǎn)云集合中的點(diǎn)云坐標(biāo)決定

(15)

式中：sz為點(diǎn)云集合中所有點(diǎn)的總數(shù)。

Riqueme et al.(2018)將結(jié)構(gòu)面假設(shè)成為不規(guī)則多邊形，并通過計(jì)算所有頂點(diǎn)之間的距離以尋求最遠(yuǎn)距離，此時(shí)的最遠(yuǎn)距離即看作該結(jié)構(gòu)面的跡長，由式(16)計(jì)算。

TL=max：ed(pipj)

=‖pi-pj‖2?i，j∈{1，…，sz}TL，ed∈，

sz∈N

(16)

同時(shí)，定義結(jié)構(gòu)面間距為同組中兩個(gè)相鄰結(jié)構(gòu)面之間的法線距離。當(dāng)結(jié)構(gòu)面理想地彼此平行時(shí)，間距可以通過式(17)～式(19)計(jì)算。

結(jié)構(gòu)面1：

Ax+By+Cz+D1=0

(17)

結(jié)構(gòu)面2：

Ax+By+Cz+D2=0

(18)

(19)

式中：結(jié)構(gòu)面1和結(jié)構(gòu)面2是兩個(gè)相鄰結(jié)構(gòu)面的平面多項(xiàng)式方程，由單位法向量(A，B，C)和常數(shù)參數(shù)D1和D2定義，SP12是這兩個(gè)結(jié)構(gòu)面的間距值。

圖6 結(jié)構(gòu)面識(shí)別結(jié)果Fig.6 Structural surface recognition resultsa.Ⅲ-B區(qū)120平臺(tái)露頭結(jié)構(gòu)面分組；b.赤平投影(等面積網(wǎng))

3.4 結(jié)構(gòu)面識(shí)別解譯結(jié)果

利用上述結(jié)構(gòu)面信息識(shí)別算法，對Ⅲ-B區(qū)域獲取到無人機(jī)點(diǎn)云數(shù)據(jù)進(jìn)行結(jié)構(gòu)面智能識(shí)別與解譯，結(jié)構(gòu)面識(shí)別結(jié)果與赤平投影圖如圖6所示。Ⅲ-B區(qū)120m平臺(tái)的巖體露頭存在4組優(yōu)勢結(jié)構(gòu)面：J1：268°∠76°，J2：310°∠40°，J3：115°∠71°，J4：344°∠78°。各組結(jié)構(gòu)面解譯后的產(chǎn)狀、跡長、間距、密度等幾何信息進(jìn)行統(tǒng)計(jì)并分析，如表 2所示。

表 2 結(jié)構(gòu)面統(tǒng)計(jì)信息Table2 Statistics of discontinuities

對于相同測線位置的結(jié)構(gòu)面識(shí)別，算法識(shí)別結(jié)果和人工量測結(jié)果是很相近的(圖7)。由于人工測量主觀性強(qiáng)，受現(xiàn)場客觀條件影響較大，很多結(jié)構(gòu)面信息無法準(zhǔn)確提取，因此機(jī)器算法所識(shí)別的結(jié)構(gòu)面數(shù)更多。表 3對比展示了兩種方法的量測結(jié)果，算法識(shí)別結(jié)果和人工量測結(jié)果傾向傾角角度差值在2°～6°之間；4組結(jié)構(gòu)面的跡長差值在0.14～0.34m之間，間距差值在0.028～0.26m之間。差值均在可接受范圍內(nèi)，證明基于點(diǎn)云數(shù)據(jù)的結(jié)構(gòu)面識(shí)別的精準(zhǔn)程度是可靠的。

表 3 算法識(shí)別與人工量測結(jié)果參數(shù)對比Table 3 Comparison of result parameters between algorithm recognition and manual measurement

圖7 兩種測量結(jié)果對比Fig.7 Comparison of the measurement results between two methodsa.邊坡測線編錄剖面；b.人工精測線量測結(jié)果；c.算法識(shí)別結(jié)果

4 三維離散元建模流程

4.1 數(shù)值模型建立

作為常用三維離散元軟件，3DEC無法直接導(dǎo)入點(diǎn)云數(shù)據(jù)進(jìn)行復(fù)雜地質(zhì)體的三維建模。本文借助Rhino軟件的三維模型構(gòu)建技術(shù)，結(jié)合Griddle軟件的網(wǎng)格處理功能實(shí)現(xiàn)任意復(fù)雜地質(zhì)體的數(shù)值網(wǎng)格模型快速剖分。整個(gè)工作過程均支持交互式界面操作，且網(wǎng)格剖分過程體現(xiàn)了高度的自動(dòng)化特點(diǎn)。因此通過“Rhino-Griddle-3DEC”聯(lián)合建模方法，可將點(diǎn)云數(shù)據(jù)導(dǎo)入3DEC完成三維數(shù)值建模。具體流程如下：

(1)確定目標(biāo)對象。本文選取Ⅲ-B區(qū)120m平臺(tái)長18m，高6m的點(diǎn)云模型開展三維數(shù)值建模(圖8a)。

圖8 建模具體流程圖Fig.8 Specific flow chart of modeling

(2)初始點(diǎn)云的稀疏處理。在3DEC中將離散裂隙網(wǎng)絡(luò)(DFN)插入數(shù)值模型時(shí)，過密的模型網(wǎng)格會(huì)導(dǎo)致DFN切割模型失敗，因此需要通過對初始點(diǎn)云進(jìn)行稀疏及優(yōu)化預(yù)處理(圖8b)，形成適合后期處理的點(diǎn)云模型。

(3)三角形網(wǎng)格生成與修復(fù)。將點(diǎn)云導(dǎo)入Geomagic Studio，點(diǎn)云中的各點(diǎn)以三角網(wǎng)格形式通過相鄰邊連接起來，達(dá)到緊密配合的效果。由于上述過程中可能會(huì)產(chǎn)生漏洞、尖銳特征和自相交三角形等錯(cuò)誤、多余連接，影響隨后的曲面生成和后續(xù)計(jì)算報(bào)錯(cuò)(圖8c)。因此在Geomagic Studio中進(jìn)行封裝后需要對曲面進(jìn)行漏洞修補(bǔ)、平滑及自相交修復(fù)，并將修復(fù)后的網(wǎng)格模型導(dǎo)入Rhino中(圖8d)。

圖9 離散網(wǎng)絡(luò)模型切割后的邊坡計(jì)算模型Fig.9 Slope calculation model after discrete network model cuttinga.初始離散裂隙網(wǎng)絡(luò)；b.簡化后離散裂隙網(wǎng)絡(luò)；c.3DEC邊坡計(jì)算模型

(4)三維實(shí)體模型生成。借助Rhino中的RhinoResurf插件完成非均勻有理B樣條(NURBS)曲面重構(gòu)及裁剪(慈瑞梅等，2004)，利用Rhino中的“擠出曲面”和“布爾運(yùn)算切割”功能，根據(jù)現(xiàn)場調(diào)查情況設(shè)置5m的厚度和底座，形成三維實(shí)體模型(圖8e)。

(5)3DEC數(shù)值模型生成。借助Griddle做水密處理，形成符合不透水定義的封閉曲面后，通過網(wǎng)格剖分得到三維數(shù)值模型，最終輸出為3DEC可以識(shí)別的GVol.3DDAT文件格式并將其導(dǎo)入3DEC軟件。如圖8f所示，3DEC中邊坡數(shù)值模型尺寸為18m8m8m，共計(jì)生成1737個(gè)離散單元體。

4.2 離散裂隙網(wǎng)絡(luò)生成

3DEC中具有強(qiáng)大的離散裂隙網(wǎng)絡(luò)(DFN)生成能力。為使DFN生成的結(jié)構(gòu)面圓盤可以成功切割邊坡模型且完成計(jì)算，離散裂隙網(wǎng)絡(luò)的范圍應(yīng)稍大于邊坡模型，尺寸為20m10m10m。DFN模型特征只受結(jié)構(gòu)面各參數(shù)統(tǒng)計(jì)分布的影響，包括裂隙尺寸(半徑)、產(chǎn)狀、位置和密度的分布統(tǒng)計(jì)(Zheng et al.，2014；Kong et al.，2020)。

表 4 Ⅲ-B區(qū)巖體結(jié)構(gòu)面參數(shù)統(tǒng)計(jì)Table 4 Statistics of rock mass discontinuities parameters in Ⅲ-B area

4.3 計(jì)算結(jié)果及分析

3DEC中的本構(gòu)模型分為塊體本構(gòu)模型和節(jié)理本構(gòu)模型，其中巖石材料設(shè)定為理想彈塑性模型，采用Mohr-Coulomb的本構(gòu)關(guān)系(Itasca，2003)。

表 5 Ⅲ-B區(qū)120￣￣m平臺(tái)邊坡巖體力學(xué)參數(shù)Table 5 Mechanical parameters of rock mass in 120m platform of the Ⅲ-B area

上部邊界為自由邊界，模型四周及底面邊界為約束邊界。初始地應(yīng)力天然工況下只考慮自重。根據(jù)巖石力學(xué)實(shí)驗(yàn)數(shù)據(jù)和現(xiàn)場工程資料，Ⅲ-B區(qū)120m平臺(tái)邊坡僅出露千枚巖，未見其他巖性出露，表 5給出了具體的巖石力學(xué)參數(shù)。

由圖10可知，模型計(jì)算到6315步時(shí)達(dá)到平衡，邊坡在天然工況下穩(wěn)定性良好。位移集中在左側(cè)陡坡邊緣處，最大位移為0.82mm。變形主要發(fā)生在Z軸(豎向)方向，表現(xiàn)為結(jié)構(gòu)面交切形成的楔形體移動(dòng)，以坡體的沉降為主。在X方向，最大位移發(fā)生在坡體臨空面，表現(xiàn)為結(jié)構(gòu)面交切形成的楔形體移動(dòng)，最大位移為0.18mm(方向沿X軸正方向)。在Y方向上，最大位移為0.3mm，為上部坡體向臨空面的滑動(dòng)。綜上可知，Ⅲ-B區(qū)120m平臺(tái)邊坡的主要位移為豎直方向上的沉降。

圖10 3DEC計(jì)算至10000步結(jié)果Fig.10 Results when 3DEC calculates to 10000 stepsa.整體位移云圖；b.X向位移云圖；c.Y向位移云圖； d.Z向位移云圖

在邊坡模型上設(shè)置6個(gè)關(guān)鍵點(diǎn)，監(jiān)測其位移隨時(shí)間的變化規(guī)律。圖11中的曲線可以看出，監(jiān)測點(diǎn)的Z向位移曲線都在相對較短的時(shí)間內(nèi)達(dá)到穩(wěn)定狀態(tài)，不隨計(jì)算時(shí)步變化，即位移時(shí)程曲線是收斂的，邊坡處于穩(wěn)定狀態(tài)。根據(jù)離散元強(qiáng)度折減法，利用3DEC的solve fos命令進(jìn)一步計(jì)算安全系數(shù)，到12363步時(shí)達(dá)到臨界破壞，此時(shí)計(jì)算安全系數(shù)所得為1.27，驗(yàn)證了天然工況下Ⅲ-B區(qū)120m平臺(tái)邊坡穩(wěn)定性狀態(tài)良好。

圖11 邊坡模型監(jiān)測點(diǎn)布置與監(jiān)測結(jié)果Fig.11 Monitoring points arrangement and the resultsa.監(jiān)測點(diǎn)位置分布；b.監(jiān)測點(diǎn)位移-時(shí)步曲線

5 結(jié) 論

以江西德興銀山銅礦露天高邊坡作為工程背景，結(jié)合無人機(jī)測量與“Rhino-Griddle-3DEC”聯(lián)合建模方法，利用結(jié)構(gòu)面統(tǒng)計(jì)信息構(gòu)建離散裂隙網(wǎng)絡(luò)切割邊坡模型，實(shí)現(xiàn)高陡巖質(zhì)邊坡穩(wěn)定性的三維離散元數(shù)值計(jì)算。主要結(jié)論如下：

(1)基于SfM三維重構(gòu)技術(shù)建立露天礦山邊坡的三維點(diǎn)云模型。通過計(jì)算點(diǎn)法向量、聚類分割等算法建立巖體結(jié)構(gòu)幾何參數(shù)識(shí)別與解譯方法。對比現(xiàn)場人工測量結(jié)果驗(yàn)證了該方法的可靠性。對優(yōu)勢分組，優(yōu)勢產(chǎn)狀，跡長分布，間距分布等結(jié)構(gòu)面信息進(jìn)行統(tǒng)計(jì)分析，構(gòu)建巖體結(jié)構(gòu)表征模型。

(2)提出“Rhino-Griddle-3DEC”的聯(lián)合建模方法，將三維點(diǎn)云模型導(dǎo)入3DEC中構(gòu)建邊坡三維離散元數(shù)值模型。根據(jù)巖體結(jié)構(gòu)表征模型建立離散裂隙網(wǎng)絡(luò)DFN，開展邊坡穩(wěn)定性計(jì)算。

(3)本文提出的結(jié)構(gòu)面識(shí)別解譯方法量測的區(qū)域大，且測量過程不影響邊坡生產(chǎn)作業(yè)。結(jié)構(gòu)面解算精度符合后期進(jìn)行巖體質(zhì)量分級或穩(wěn)定性計(jì)算的要求，顯著提高了大規(guī)模高邊坡的巖體編錄效率。結(jié)合邊坡三維點(diǎn)云建模與巖體結(jié)構(gòu)面統(tǒng)計(jì)分析，“Rhino-Griddle-3DEC”建模方法能夠快速建立巖體結(jié)構(gòu)的三維離散元模型，獲得更加符合實(shí)際的數(shù)值模擬計(jì)算結(jié)果。