巖體結(jié)構(gòu)面面密度的數(shù)字化統(tǒng)計方法及其應(yīng)用

鄭 健，章楊松，李曉昭，石杏喜

(1.南京理工大學(xué)土木工程系，江蘇南京210094;2.南京大學(xué)地球科學(xué)與工程學(xué)院，江蘇南京210093)

巖體結(jié)構(gòu)面是控制巖體強弱和穩(wěn)定的主要因素，在工程和理論研究中對巖體結(jié)構(gòu)面進行統(tǒng)計和分析尤為重要。巖體結(jié)構(gòu)面密度可表征結(jié)構(gòu)面分布特征，是其三維網(wǎng)絡(luò)模擬的基本參數(shù)。結(jié)構(gòu)面密度在不同空間維度上，有不同定義方式，一般分為線密度、面密度、體密度［1］，其間在一定條件下可相互轉(zhuǎn)換［2］。其中，面密度定義為單位面積巖體內(nèi)所包含結(jié)構(gòu)面跡線中點數(shù)，或定義為單位面積內(nèi)所有結(jié)構(gòu)面長度之和［3］，本文面密度采用第一種定義。

傳統(tǒng)的獲取結(jié)構(gòu)面信息的方法包括測線法及測窗法，其測量及記錄依靠人工進行。近年來，應(yīng)用三維激光掃描技術(shù)及近景攝影測量技術(shù)獲取結(jié)構(gòu)面信息的研究也得到了發(fā)展［4～7］。其中，應(yīng)用三維激光掃描技術(shù)的后處理軟件也陸續(xù)得以開發(fā)，但還存在一些問題，如掃描點云數(shù)據(jù)量巨大、結(jié)構(gòu)面識別半自動化等;數(shù)字近景攝影測量技術(shù)中結(jié)構(gòu)面信息的解譯及統(tǒng)計系統(tǒng)也還不夠成熟。王鵬等［8］利用GPS-RTK技術(shù)量測甘肅北山裂隙信息，并以ArcGIS平臺建立裂隙屬性數(shù)據(jù)庫，采用地質(zhì)統(tǒng)計分析方法探索裂隙空間分布特征。這種方法可得到裂隙水平投影面的面密度值及水平面內(nèi)的變化規(guī)律。

本文利用GPS-RTK技術(shù)詳細(xì)調(diào)查記錄了甘肅北山研究區(qū)域大量花崗巖巖體露頭情況，獲得了豐富的裂隙跡線測量記錄。以此為基礎(chǔ)數(shù)據(jù)，開發(fā)了一種自動統(tǒng)計估算結(jié)構(gòu)面面密度的數(shù)字化后處理方法。該方法充分利用地形信息，可得到具體產(chǎn)狀巖體露頭面(或測窗)結(jié)構(gòu)面面密度值。以兩個測區(qū)數(shù)據(jù)為例，闡述此方法并說明方法的高效性與應(yīng)用的多樣性。

1 測區(qū)巖體結(jié)構(gòu)面發(fā)育基本特征

甘肅北山地區(qū)是我國高放廢物地質(zhì)處置庫建造的預(yù)選場址之一。高放廢物地質(zhì)處置場址的巖體結(jié)構(gòu)面發(fā)育特征的系統(tǒng)研究尤為重要。巖體結(jié)構(gòu)面的調(diào)查研究內(nèi)容主要包括結(jié)構(gòu)面的空間分布、密度、連通程度等特征［9］。在此，重點闡述巖體結(jié)構(gòu)面面密度方面的研究成果。

測區(qū)Ⅰ位于甘肅北山地區(qū)舊井地段十月井?dāng)嗔褞鞅狈较?，包?個露頭面，共獲取4 224條有效裂隙跡線，調(diào)查區(qū)域面積約為7 300 m2;測區(qū)Ⅱ位于芨芨槽塊段，包含2個露頭面，共獲取1 075條有效裂隙跡線，調(diào)查區(qū)域面積約為5 300 m2。根據(jù)人工測量產(chǎn)狀進行分組，各區(qū)分組結(jié)果見表1。

表1 測區(qū)結(jié)構(gòu)面分組情況Table 1 Dominant partitioning of discontinuities

2 數(shù)字化統(tǒng)計方法

2.1 建立跡線三維模型

跡線三維模型是進行研究的基礎(chǔ)，所建模型相當(dāng)于真實露頭的數(shù)字化樣本，在此基礎(chǔ)上的各項研究具有可重復(fù)性。同時，包含結(jié)構(gòu)面信息的數(shù)字化模型，可以借助計算機來完成以往人工統(tǒng)計工作，因此，大大提高了結(jié)構(gòu)面信息的處理效率和準(zhǔn)確性。

模型建立的基礎(chǔ)數(shù)據(jù)為露頭面跡線控制點(端點或明顯轉(zhuǎn)折點)坐標(biāo)，可以GPS-RTK技術(shù)獲得。在實地調(diào)查中，跡線以折線(或直線)形式測量，而在建模及統(tǒng)計過程中，為求簡便，以首尾端點所連直線段表示。本模型及其后研究對地形精度要求不高，可根據(jù)一定數(shù)量的跡線控制點形成數(shù)字地面。以研究區(qū)所有裂隙跡線控制點做為散亂點，先形成 TIN模型(Triangulated Irregular Network，不規(guī)則三角網(wǎng))，再以網(wǎng)格面擬合成光滑曲面形成地面，將跡線(空間直線)與數(shù)字地面結(jié)合，便得到了跡線三維模型。圖1即為測區(qū)Ⅰ、Ⅱ跡線三維模型。測區(qū)Ⅰ中的4個露頭面分別記為A，B(由B1，B2組成)，C，D 露頭;測區(qū)Ⅱ中的2個露頭面分別記為E，F(xiàn)露頭。圖1中坐標(biāo)系為空間直角坐標(biāo)系，y軸正方向為北方向，x軸正方向為東方向，z軸為鉛垂軸;圖中網(wǎng)格面為擬合數(shù)字地面，不同顏色的直線線條代表不同組裂隙跡線;地面上無跡線區(qū)域為未調(diào)查區(qū)域。

2.2 數(shù)字化測窗的布置與原理

建立跡線三維模型后，可由自編程序在研究區(qū)域自動布置矩形測窗，以測區(qū)ⅠA、B露頭面區(qū)域跡線三維模型為例(圖2)。圖2中，同一顏色網(wǎng)格面代表一個矩形測窗，未調(diào)查區(qū)域不布置測窗，顯示為灰色網(wǎng)格地面。

這些地面矩形測窗并非現(xiàn)有估算理論中的二維平面測窗(圖3)，說明地面矩形測窗與程序計算測窗之間的關(guān)系。程序根據(jù)地面矩形測窗ABCD的4角點坐標(biāo)，以最小二乘原理擬合出一平面，在此平面內(nèi)的測窗(即A'B'C'D')即為估算理論中的矩形測窗。擬合平面上矩形測窗4個角點的x，y值與地面矩形測窗的4個角點x，y值相同。

下文中所提矩形、圓形測窗，若無特殊說明，均為擬合平面上的測窗。

圖1 測區(qū)Ⅰ、Ⅱ跡線三維模型Fig.1 3D digital traces model of regionⅠandⅡ

圖2 測區(qū)ⅠA、B露頭地面測窗布置情況Fig.2 Disposition of sampling windows on the surface

圖3 矩形、圓形測窗與地面測窗關(guān)系示意圖Fig.3 Relationship between sampling windows on the surface and windows used in the program

為了更好地進行圓形測窗與矩形測窗面密度估算結(jié)果的對比，各圓形測窗布置在矩形測窗所在擬合平面上，圓心O取為矩形測窗對角線交點，圓形測窗面積與矩形測窗面積平均值相同。實際上，在程序完成地面矩形測窗的自動布置后，擬合平面上矩形、圓形測窗的大小、位置信息可通過程序計算得到，在圖2中未顯示出來。

2.3 面密度統(tǒng)計估算程序算法

2.3.1 估算理論的選取

測窗布置完成后，需選取估算方法對各個測窗內(nèi)的面密度進行估算。目前國際上對于面密度估算的方法主要有2種:一種是Kulatilake等的面密度估算方法［10］，另一種是 Mauldon 估算方法［11］。Mauldon 估算面密度方法基于凸多邊形測窗，圓形、矩形測窗都適用。相比Mauldon方法，Kulatilake等估算方法計算更復(fù)雜，同時需要獲取研究區(qū)域跡長的概率密度函數(shù)。程序采用Mauldon估算方法，同時應(yīng)用矩形及圓形測窗，進行研究區(qū)域面密度的估算。

Mauldon提出的結(jié)構(gòu)面面密度估算以下式計算:

式中:λ——結(jié)構(gòu)面面密度;

A——測窗面積;

N0，N1，N2——窗口中兩端均不可見、一端可見、兩端可見類型的跡線條數(shù)(同時，為表述方便，也作為各類型表示符號);N=N0+N1+N2，代表以上三種跡線數(shù)目總和。

由式(1)可知，此估算方法的缺點是忽略了貫穿型跡線對面密度的影響，這一缺點可通過適當(dāng)增大測窗面積以減少N0來彌補。

2.3.2 程序算法

由式(1)可知，估算值取決于測窗面積及N1，N2，其中，測窗面積可由程序自動計算得到。在三維空間中N1，N2的計數(shù)較困難，可將所有元素正投影到水平面上(z=0平面)。由于正投影不改變跡線與測窗的交切關(guān)系，所以在投影面上做跡線類型判斷不改變原有結(jié)果，三維問題由此轉(zhuǎn)化為平面不同類型跡線的統(tǒng)計問題。

由于跡線兩端點坐標(biāo)均已知，各擬合平面產(chǎn)狀及平面上矩形、圓形測窗的尺寸及位置信息可由程序自動計算得到，根據(jù)跡線兩端點投影坐標(biāo)與測窗投影圖形位置的相互關(guān)系可判斷跡線類型。由此統(tǒng)計得各測窗N0，N1，N2計數(shù)，借由式(1)可求得各測窗面密度估算值。同時，由于本研究中跡線全跡長已知，可借投影方法判斷跡線中點與相應(yīng)測窗的關(guān)系，由程序計算得到測窗的真實面密度值。

3 方法應(yīng)用結(jié)果及分析

3.1 面密度統(tǒng)計估算結(jié)果對比

將上述數(shù)字化統(tǒng)計方法應(yīng)用于甘肅北山13個測區(qū)，得到了2 000余測窗的統(tǒng)計結(jié)果。在此，以測區(qū)Ⅰ、Ⅱ的統(tǒng)計結(jié)果為例，說明應(yīng)用情況。

以各測窗面密度真值λt與估算值λre(矩形測窗面密度估算值)、λce(圓形測窗面密度估算值)對比。圖4為測區(qū)Ⅰ(共189個測窗)第1組結(jié)構(gòu)面面密度對比結(jié)果。測區(qū)Ⅰ剩余組結(jié)構(gòu)面、測區(qū)Ⅱ的各組結(jié)構(gòu)面面密度對比圖與圖4類似，不再列出。由圖4可觀察到λre，λce雖然在個別窗口與真值λt相差較大，但就研究區(qū)域面密度的變化情況來看，兩者變化趨勢與面密度真值變化趨勢是十分一致的。

為評價矩形測窗及圓形測窗的估算精度，需計算估算誤差:

式中:λe——面密度估算值。

圖4 測區(qū)Ⅰ第1組結(jié)構(gòu)面面密度估算結(jié)果與真值對比Fig.4 Comparison between the estimated and true values of 2D density of Set 1 in regionⅠ

應(yīng)用式(2)時將λre，λce值代入計算即可。誤差樣本均值記為 eλr，樣本標(biāo)準(zhǔn)差記為 Sλr。圓形測窗估算誤差樣本均值記為 eλc，樣本標(biāo)準(zhǔn)差記為Sλc。計算結(jié)果見表2。

由表2可知，在測窗尺寸合適的情況下，Mauldon估算理論的估算精度在總體上是能滿足工程需求的。應(yīng)用矩形或圓形測窗進行估算，其誤差樣本均值及標(biāo)準(zhǔn)差均無較大差異，說明在Mauldon估算理論下，應(yīng)用矩形或圓形測窗對估算結(jié)果影響不大。但在應(yīng)用矩形測窗時，相比圓形測窗，連續(xù)布置的矩形測窗可無間隙、無重復(fù)的覆蓋研究區(qū)域，更有利于掌握面密度的地域分布情況。

表2 測窗估算誤差樣本的均值及標(biāo)準(zhǔn)差Table 2 Mean and standard deviation values of the estimating error

3.2 研究區(qū)面密度分布情況

描述結(jié)構(gòu)面的幾何特征參數(shù)如間距、跡長等，其調(diào)查、研究成果豐富，參數(shù)的概率統(tǒng)計模型一般服從負(fù)指數(shù)分布、對數(shù)正態(tài)分布或正態(tài)分布。但是，鮮有對結(jié)構(gòu)面面密度統(tǒng)計模型分布的研究，其主要原因是難以獲取大范圍、充足數(shù)量的實測面密度數(shù)據(jù)。同時，若對某一區(qū)域進行面密度分布研究，理想情況下的測窗布置應(yīng)是連續(xù)、無重復(fù)的覆蓋研究區(qū)域。本文數(shù)字化統(tǒng)計方法中，程序可自動布置矩形測窗，不僅滿足連續(xù)無重復(fù)布置要求，還可任意調(diào)整測窗大小。

以各測窗真實面密度λt為對象進行統(tǒng)計。由于程序計算得到的面密度其所屬的巖體露頭(或測窗)產(chǎn)狀各不相同，需向特定產(chǎn)狀剖面統(tǒng)一［12］，本文將面密度向與優(yōu)勢產(chǎn)狀結(jié)構(gòu)面垂直的剖面統(tǒng)一:

式中:λs——統(tǒng)一產(chǎn)狀巖體剖面上結(jié)構(gòu)面面密度;

δs——統(tǒng)一產(chǎn)狀巖體剖面與優(yōu)勢產(chǎn)狀結(jié)構(gòu)面夾角，文中取為90°;

λ0——測窗真實結(jié)構(gòu)面面密度;

δ0——該測窗與本組優(yōu)勢產(chǎn)狀結(jié)構(gòu)面夾角。

完成各區(qū)各組結(jié)構(gòu)面面密度的修正后，進行直方圖統(tǒng)計，再觀察其分布規(guī)律，以合適的函數(shù)擬合其分布，如圖5。圖5為兩種代表性分布形式:圖5(a)為測區(qū)Ⅰ第1組結(jié)構(gòu)面面密度分布，以對數(shù)正態(tài)函數(shù)擬合合適，此分布在測區(qū)Ⅰ，Ⅱ中僅此1例;圖5(b)為測區(qū)Ⅱ第1組結(jié)構(gòu)面面密度分布，以雙參數(shù)負(fù)指數(shù)分布函數(shù)［13］擬合合適，測區(qū)Ⅰ第2，3組結(jié)構(gòu)面及測區(qū)Ⅱ的第2，3，4組結(jié)構(gòu)面面密度分布形式與其相同。

雙參數(shù)負(fù)指數(shù)分布具體形式如下:

式中:a，b——擬合曲線待定參數(shù)。

圖5 結(jié)構(gòu)面面密度頻率統(tǒng)計直方圖與擬合曲線及公式Fig.5 Frequency histogram，fitting curve and formula of discontinuities

各區(qū)各分組面密度參數(shù)及擬合分布參數(shù)見表3，R2為在0～1之間的判定系數(shù)，越接近1說明擬合效果越好。測區(qū)Ⅰ第1組擬合參數(shù)在圖5(a)中。

表3 測區(qū)面密度分布情況及擬合參數(shù)Table 3 Distribution information of 2D density and fitting parameters

由表3可觀察到，測區(qū)Ⅰ各組面密度均值遠(yuǎn)大于測區(qū)Ⅱ，這是由于測區(qū)Ⅰ臨近十月井?dāng)嗔褞?，總體上裂隙更為發(fā)育。實際上，根據(jù)統(tǒng)計結(jié)果，甘肅北山芨芨槽塊段另11個測區(qū)(非臨近斷裂帶區(qū)域)面密度分布情況均與測區(qū)Ⅱ相近，以雙參數(shù)負(fù)指數(shù)分布函數(shù)可以取得良好的擬合效果，R2均大于0.9。

4 結(jié)論

本文所提出的結(jié)構(gòu)面面密度的數(shù)字化統(tǒng)計方法是一種結(jié)構(gòu)面信息后處理方法，應(yīng)用該方法的基礎(chǔ)資料為裂隙跡線相對坐標(biāo)，文中測區(qū)Ⅰ，Ⅱ數(shù)據(jù)以GPSRTK技術(shù)獲得。在建立測區(qū)跡線三維模型的基礎(chǔ)上，所作的主要工作及結(jié)論有:

(1)實現(xiàn)了數(shù)字化測窗的連續(xù)自動化布置，編制了統(tǒng)計估算結(jié)構(gòu)面面密度的程序。

(2)獲取了測區(qū)Ⅰ，Ⅱ面密度估算值與真值，通過對比，驗證了Mauldon面密度估算理論的準(zhǔn)確性。

(3)獲取了測區(qū)Ⅰ，Ⅱ共計291個測窗面密度數(shù)據(jù)，以此研究其分布規(guī)律。結(jié)合實際13個測區(qū)統(tǒng)計結(jié)果，發(fā)現(xiàn)絕大多數(shù)測區(qū)其結(jié)構(gòu)面面密度分布以雙參數(shù)負(fù)指數(shù)分布函數(shù)擬合的效果良好。

該方法在處理大范圍、數(shù)據(jù)量巨大的巖體結(jié)構(gòu)面面密度信息的工作中具有較大優(yōu)勢，未來，還可以以跡線三維模型為基礎(chǔ)，進行結(jié)構(gòu)面間距計算、估算連通率、進行面密度空間變異性等研究，應(yīng)用前景廣闊。

［1］ 王輝，錢海濤.節(jié)理化巖體結(jié)構(gòu)面分布密度的確定方法及其應(yīng)用［J］.地球與環(huán)境，2005，33(增刊1):124-128.［WANG H，QIAN H T.Thetechnique application ofassessing thedistribution densityof structures in rock mass［J］.Earth and Environment，2005，33(Sup 1):124-128.(in Chinese)］

［2］ 伍法權(quán).統(tǒng)計巖體力學(xué)原理［M］.武漢:中國地質(zhì)大學(xué)出版社，1993.［WU F Q.Principles of statistical mechanics of rock masses［M］.Wuhan:China University of Geosciences Press，1993.(in Chinese)］

［3］ 趙文，林韻梅.結(jié)構(gòu)面密度新指標(biāo)［J］.中國礦業(yè)，1994，3(3):42-44.［ZHAO W，LIN Y M.A new parameter of discontinuity density［J］.China Mining Magazine，1994，3(3):42-44.(in Chinese)］

［4］ 霍俊杰，黃潤秋，董秀軍，等.3D激光掃描與巖體結(jié)構(gòu)精細(xì)測量方法比較研究［J］.湖南科技大學(xué)學(xué)報(自然科學(xué)版)，2011，26(3):39-44.［HUO J J，HUANG R Q，DONG X J，et al.Comparison analysis on 3D laser scanning technology and refined mesh investigation of rock mass structure［J］.Journal of Hunan University of Science＆Technology(Natural Science Edition)，2011，26(3):39-44.(in Chinese)］

［5］ 劉昌軍，高立東，丁留謙，等.應(yīng)用激光掃描技術(shù)進行巖體結(jié)構(gòu)面的半自動統(tǒng)計研究［J］.水文地質(zhì)工程地質(zhì)，2011，38(2):52-57.［LIU C J，GAO L D，DING L Q，et al.Research of semi-automatic statistics of rock mass discontinuity applying laser scanning technology［J］.Hydrogeology ＆ Engineering Geology，2011，38(2):52-57.(in Chinese)］

［6］ 王鳳艷.數(shù)字近景攝影測量快速獲取巖體裂隙信息的工程應(yīng)用［D］.長春:吉林大學(xué)，2006.［WANG F Y.Engineeringapplication ofrapid acquiring rock mass fractures information with digital close range photogrammetry［D］.Changchun:Jilin University，2006.(in Chinese)］

［7］ 王鳳艷，陳劍平，楊國東，等.基于數(shù)字近景攝影測量的巖體結(jié)構(gòu)面幾何信息結(jié)算模型［J］.吉林大學(xué)學(xué)報(地球科學(xué)版)，2012，42(6):1839-1846.［WANG F Y，CHEN J P，YANG G D，et al.Solution model of Geometrical information of rock mass discontinuity based on digital close rang photogrammetry［J］.Journal of Jilin University(Earth Science Edition)，2012，42(6):1839-1846.(in Chinese)］

［8］ 王鵬，李曉昭，章楊松，等.基于GIS的甘肅北山花崗巖裂隙密度地質(zhì)統(tǒng)計分析［J］.工程地質(zhì)學(xué)報，2013，21(1):115-122.［WANG P，LI X Z，ZHANG Y S，et al.GIS based geostatistical analysis of fracture density of granite rock in Beishan area Gansu Province［J］.Journal of Engineering Geology，2013，21(1):115-122.(in Chinese)］

［9］ 郭永海，王駒.高放廢物地質(zhì)處置中的地質(zhì)、水文地質(zhì)、地球化學(xué)關(guān)鍵科學(xué)問題［J］.巖石力學(xué)與工程學(xué)報，2007，26(增刊 2)，3926-3931.［GUO Y H，WANG J.Key scientific issues of geology，hydrogeology and geochemistry in high-level radioactive waste geological disposal［J］.Chinese Journal of Rock Mechanics and Engineering，2007，26(Sup 2)，3926-3931.(in Chinese)］

［10］ Kulatilake P H S W，Wu T H.The density of discontinuity traces in sampling windows［J］.International Journal of Rock Mechanics and Mining Sciences，1984，21:345-347.

［11］ Mauldon M.Estimating mean fracture trace length and density from observations in convex windows［J］.Rock Mechanics and Rock Engineering，1998，31(4):201-216.

［12］ 金曲生，范建軍，王思敬.結(jié)構(gòu)面密度計算法及其應(yīng)用［J］.巖石力學(xué)與工程學(xué)報，1998，17(3):273-278.［JING Q S，F(xiàn)AN J J，WANG S J.A new method to calculate density of discontinuity and its application［J］.Chinese Journal of Rock Mechanics and Engineering，1998，17(3):273-278.(in Chinese)］

［13］ 胡秀宏，伍法權(quán).巖體結(jié)構(gòu)面間距的雙參數(shù)負(fù)指數(shù)分布研究［J］.巖土力學(xué)，2009，30(8):2353-2358.［HU X H，WU F Q.Research on two-parameter negative exponential distribution of discontinuity spacing in rock mass［J］.Rock and Soil Mechanics，2009，30(8):2353-2358.(in Chinese)］