巖體結構面面密度的數字化統計方法及其應用

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

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

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

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

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

1 測區巖體結構面發育基本特征

甘肅北山地區是我國高放廢物地質處置庫建造的預選場址之一。高放廢物地質處置場址的巖體結構面發育特征的系統研究尤為重要。巖體結構面的調查研究內容主要包括結構面的空間分布、密度、連通程度等特征［9］。在此，重點闡述巖體結構面面密度方面的研究成果。

測區Ⅰ位于甘肅北山地區舊井地段十月井斷裂帶西北方向，包含4個露頭面，共獲取4 224條有效裂隙跡線，調查區域面積約為7 300 m2;測區Ⅱ位于芨芨槽塊段，包含2個露頭面，共獲取1 075條有效裂隙跡線，調查區域面積約為5 300 m2。根據人工測量產狀進行分組，各區分組結果見表1。

表1 測區結構面分組情況Table 1 Dominant partitioning of discontinuities

2 數字化統計方法

2.1 建立跡線三維模型

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

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

2.2 數字化測窗的布置與原理

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

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

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

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

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

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

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

2.3 面密度統計估算程序算法

2.3.1 估算理論的選取

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

Mauldon提出的結構面面密度估算以下式計算:

式中:λ——結構面面密度;

A——測窗面積;

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

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

2.3.2 程序算法

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

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

3 方法應用結果及分析

3.1 面密度統計估算結果對比

將上述數字化統計方法應用于甘肅北山13個測區，得到了2 000余測窗的統計結果。在此，以測區Ⅰ、Ⅱ的統計結果為例，說明應用情況。

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

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

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

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

應用式(2)時將λre，λce值代入計算即可。誤差樣本均值記為 eλr，樣本標準差記為 Sλr。圓形測窗估算誤差樣本均值記為 eλc，樣本標準差記為Sλc。計算結果見表2。

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

表2 測窗估算誤差樣本的均值及標準差Table 2 Mean and standard deviation values of the estimating error

3.2 研究區面密度分布情況

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

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

式中:λs——統一產狀巖體剖面上結構面面密度;

δs——統一產狀巖體剖面與優勢產狀結構面夾角，文中取為90°;

λ0——測窗真實結構面面密度;

δ0——該測窗與本組優勢產狀結構面夾角。

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

雙參數負指數分布具體形式如下:

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

圖5 結構面面密度頻率統計直方圖與擬合曲線及公式Fig.5 Frequency histogram，fitting curve and formula of discontinuities

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

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

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

4 結論

本文所提出的結構面面密度的數字化統計方法是一種結構面信息后處理方法，應用該方法的基礎資料為裂隙跡線相對坐標，文中測區Ⅰ，Ⅱ數據以GPSRTK技術獲得。在建立測區跡線三維模型的基礎上，所作的主要工作及結論有:

(1)實現了數字化測窗的連續自動化布置，編制了統計估算結構面面密度的程序。

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

(3)獲取了測區Ⅰ，Ⅱ共計291個測窗面密度數據，以此研究其分布規律。結合實際13個測區統計結果，發現絕大多數測區其結構面面密度分布以雙參數負指數分布函數擬合的效果良好。

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

［1］ 王輝，錢海濤.節理化巖體結構面分布密度的確定方法及其應用［J］.地球與環境，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］ 伍法權.統計巖體力學原理［M］.武漢:中國地質大學出版社，1993.［WU F Q.Principles of statistical mechanics of rock masses［M］.Wuhan:China University of Geosciences Press，1993.(in Chinese)］

［3］ 趙文，林韻梅.結構面密度新指標［J］.中國礦業，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激光掃描與巖體結構精細測量方法比較研究［J］.湖南科技大學學報(自然科學版)，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］ 劉昌軍，高立東，丁留謙，等.應用激光掃描技術進行巖體結構面的半自動統計研究［J］.水文地質工程地質，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］ 王鳳艷.數字近景攝影測量快速獲取巖體裂隙信息的工程應用［D］.長春:吉林大學，2006.［WANG F Y.Engineeringapplication ofrapid acquiring rock mass fractures information with digital close range photogrammetry［D］.Changchun:Jilin University，2006.(in Chinese)］

［7］ 王鳳艷，陳劍平，楊國東，等.基于數字近景攝影測量的巖體結構面幾何信息結算模型［J］.吉林大學學報(地球科學版)，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的甘肅北山花崗巖裂隙密度地質統計分析［J］.工程地質學報，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］ 郭永海，王駒.高放廢物地質處置中的地質、水文地質、地球化學關鍵科學問題［J］.巖石力學與工程學報，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］ 金曲生，范建軍，王思敬.結構面密度計算法及其應用［J］.巖石力學與工程學報，1998，17(3):273-278.［JING Q S，FAN 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］ 胡秀宏，伍法權.巖體結構面間距的雙參數負指數分布研究［J］.巖土力學，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)］