岩体结构面面密度的数字化统计方法及其应用

郑 健，章杨松，李晓昭，石杏喜

(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)］