节理产状极点分形维计算方法的程序实现

基金项目:地质灾害防治与地质环境保护国家重点实验室开放基金（SKLGP2012K032）

作者简介:宋丽娟(1983)，女，河南郑州人，博士研究生，主要从事岩土体稳定性方面的研究工作。Email：sljxinfeiyang@163.com

通讯作者:许模(1963)，男，重庆人，教授，博士生导师主要从事工程水文地质、环境地质方面的研究。Email：xm@cdut.edu.cn。DOI：10.13476/j.cnki.nsbdqk.2014.04.033

摘要：针对节理产状极点分布分形维D计算方法中的赤道圆网格剖分法，进行了程序设计。设计过程中利用节理产状极点的极坐标与产状倾向、倾角间的转化关系，合理地将节理产状极点图数值化，使用Visual Basic for Applications（VBA）进行编程，并在Excel中直接进行数据处理，简化编程，使该程序简单易用，有利于节理产状极点分形维D计算方法的推广和应用。最后，以大岗山水电站坝区左岸6组优势裂隙资料为例，计算建基岩体开挖前后各组节理产状极点分布的分形维D，与已有文献所提节理产状极点分形维D的相关结论进行对比，验证了程序的正确性与可用性。

关键词：赤平投影；不稳定块体；极坐标；分形维D

中图分类号：TU45文献标志码：A文章编号：16721683（2014）04015105

Program Implementation of Fractal Dimension Algorithm of Orientation Pole Distribution for Joints

SONG Lijuan1，XU Mo1，LU Shuqiang2

(1.State Key Laboratory of Geohazard Prevention and Geoenvironment Protection，Chengdu University of Technology，

Chengdu 610059，China；2.Key Laboratory of Geological Hazards on Three Gorges Reservoir

Area of Ministry of Education，China Three Gorges University，Yichang 443002，China)

Abstract:The distribution characteristics of joints with dominant orientation are difficult to determine because of the randomness and nonuniformity of spatial distribution of joints,and fractal dimension D is an important approach to describe quantitatively the orientation distribution using the boundary combination based on the judgment of instable block with stereographic projection.According to the mesh generation method of equatorial circle in fractal dimension D,the corresponding program was designed in this paper.In the design process,dip direction and dip angle were converted to polar coordinates and to realize the numerical representation of orientation pole distribution for joints.The program was designed with Visual Basic for Applications (VBA) and the data were processed directly in Excel,which made the program easy to operate and was beneficial to promote the application in the fractal dimension algorithm of orientation pole distribution for joints.Finally,based on the data of six sets of dominant joints collecting from the left dam foundation in the Dagangshan Hydropower Station,the fractal dimension D of six sets of dominant joints of the foundation rock was calculated before and after the excavation,which was verified with the conclusions on the fractal dimension D of orientation pole distribution for joints in previous studies,indicatingthe correctness and feasibility of the program.

Key words:stereographic projection;instable block;polar coordinate;fractal dimension D

岩体稳定性分析是目前工程地质界的一个重大难题，特别是确定不稳定块体边界面及其组合模式。目前通常利用赤平投影来确不稳定块体的边界面组成 ［16］。该方法中通常以结构面的平均产状代表结构面，对于构造断层、断层、岩脉等Ⅰ至Ⅳ级结构面来说，在某一空间范围内其几何特征相对稳定，具有一定的确定性，该方法的可靠性较高；但是，对于节理裂隙、微裂隙等Ⅴ级结构面来说，在某一空间范围内，其发育的随机性和不均匀性将导致优势节理组的产状具有一定的分散性［7］，此时以平均产状代表其节理的发育状况，其可信度较低，则该方法的可靠度较低。鉴于此，在利用赤平投影对不稳定块体的边界面进行确定的同时，需要对节理产状的分布情况进行详细研究，进而判断节理组作为不稳定块体边界组成的可靠度。

历来，国内外学者非常重视节理产状的研究，如：早期，R.Fisher［8］建立结构面产状的Fisher分布模型；C.Bingham［9］提出了关于产状分布的Bingham模型；R.E.Goodman等［1011］指出在坚硬岩体中，结构面产状服从半球正态分布；之后，一些学者尝试动态聚类分析法、k均值聚类分析法、模糊等价聚类方法、模糊C均值聚类方法、概率模型方法等［1215］对节理产状进行了深入研究；然而，上述研究成果对节理产状的描述依旧停留在定性阶段。分形理论建立之后，陈剑平等［7,16］试图利用分形维D来定量描述产状的分布特征，取得了较理想的结果。

本文拟根据节理产状极点分布分形维D的计算方法，进行程序设计，且通过引入产状极点坐标与倾向、倾角之间的转化关系，从而将节理产状极点图数值化，简化编程设计，使节理产状极点分布分形维D的计算更容易推广应用。

1原理

1.1Schmidt图等面积方格网剖分原理

目前，对Schmidt图等面积网格剖分主要有两种方法：方法1是借鉴S.M.Miller建立下半球Schmidt等面积投影的方法，通过计算机编程动态地生成下半球Schmidt等面积投影网格［7］；方法2是根据赤道圆与下半球球面呈垂直投影的关系，直接对赤道圆采用特定的网格剖分，以实现Schmidt图等面积网格剖分［16］。两种方法的原理及参数对比见表1，其原理具体参见文献［7］和文献［16］，本文主要针对方法2［16］的原理进行程序设计。

表1Schmidt图等面积方格网剖分方法对比

Table 1Comparison of two Schmidt

subdivision methods with equal area

对比项方法1方法2原理下半球Schmidt等面积投影原理对赤道圆按特定方法进行网格剖分编程设计的自变量bi、θinD的计算结果一致一致编程难易复杂简单1.2分形维D计算原理

依据方法2［16］所述原理，Schmidt圆的网格剖分方法如下见图1：沿半径R方向n等分得n个圆环，各圆环宽度均为d=R/n，中心圆为第1环，环半径等于环宽度，即r1=d1=R/n。以中心圆面积为A1=πr21=π（R/n）2基础面积，分别将各个圆环划分为面积为π（R/n）2的小方格见图2。由表2可知，当i = 1，2，3，…，n时，bi= 1，3，5，…，n2-（n-1）2，即bi为公差为2的等差数列，则n个圆环中的方格总数为NC=∑bi=n（1+bn）/2。若假定Schmidt圆面积为1个单位，该单位可取任何值，参考文献［7］中对其取值1.0×105，则小方格边长(相对)为：

δ=1.0×105NC(1)

图1等分为n个圆环的Schmidt圆

Fig.1Schmidt circle divided into n rings

表2Schmidt圆等面积方格剖分

Table 2Net subdivision with equal area of Schmidt circle

项目公式第i环环宽度di=R/n第i环环半径ri=i（R/n）第i环环面积Ai=［i2-（i-1）2］π（R/n）2第i环方格数bi=Ai/Ai=i2-（i-1）2第i环每一小方格所占角度α=360°/bi上述Schmidt圆等面积方格网的剖分过程中，仅有一个自变量n。当改变n时，相应NC改变，随之小方格边长δ改变，被节理节点占有的方格数N（δ）也改变，则可求出产状极点分形维D为：

D=lnN（δ）lnδ(2)

分形维D为lnδ-lnN（δ）曲线中无标度区间内近似直线段的曲线段斜率。

图2随n值变化的Schmidt等面积投影网

Fig.2Variation of Schmidt net of equal area with value of n

2程序实现

2.1编程原理

众所周知，节理产状极点在极点图的位置主要根据产状倾向和倾角来确定，图3给出了产状倾向320°、倾角20°的极点制作过程：首先找出倾向320°的位置，见图3中射线然后确定倾角20°所在圆位置，则射线与圆的交点即为所求产状极点，见图3。

圆中的任一点均可用一个极坐标（r,θ）来表示，则某一产状倾向、倾角与极坐标之间的关系如下：

图3产状极点确定示意图

Fig.3Schematic of orientation pole

r=R-90°-倾角90°θ=倾向(3)

根据公式（3）可以将节理产状转化为极坐标形式。

根据上述Schmidt圆等面积方格网的剖分原理和分形维D的计算原理，将各个方格极坐标域化，如第一环内的小方格域化后为{(r,θ)0≤r<2d,0°≤θ<360°}；然后根据数值化的坐标来判断各个产状极点是否包含在某一极坐标域中，求出极点占有的极坐标域的个数，即为小方格个数N(δ)。

2.2编程界面及主代码

通过极坐标转化，直接将节理产状极点图中的极点进行数据化处理，即将几何图像转为数据处理，使得编程及所用软件更加方面可靠。本文采用Visual Basic for Applications（VBA）进行编程，而在Excel中直接进行数据处理。图4是用VBA所编产状极点图分形维计算界面，其中有三个操作按钮：“Clear Worksheet”、“Import File”和“Counting”。“Clear Worksheet”用来清除Excel中的所有数据；“Import File”用来导入产状倾向、倾角数据，分别存在A列和B列；“Counting”用来计算被极点占用的小方格数，输出在result中，小方格数的占用主要与输入参数n有关，n为上述章节中的环数。右侧用来记录环数n和被极点占据的小方格数。然后利用公式（1）求出小方格边长δ，进而求出lnδ和lnN(δ)。最后将数据导入SPSS成图，并求其分形维D。

图4节理产状极点分布分形维编程设计界面

Fig.4Program interface of fractal dimension of

orientation pole distribution for joints

以上程序主代码如下：

Private Sub CommandButton3_Click()

Dim C1() As Double,C2() As Double,Z1() As Double,Z2() As Double,i As Long,j As Long,a As Long,r As Integer,n As Integer,d As Double

Dim ng As Integer,count() As Integer,mid As Integer,an As Integer

For i=1 To 100000

If (Cells(i,1) <> "") Then

a=i

Else

Exit For

End If

Next i

If Cells(a,2)="" Then

Cells(a,1)=""

a=a1

End If

ReDim C1(a),C2(a),Z1(a),Z2(a) As Double

For i=1 To a

C1(i)=Cells(i,1)

Z1(i)=C1(i)

C2(i)=Cells(i,2)

Z2(i)=90  C2(i)

Next i

n=Cells(2,3)

d=90/n

ng=n^2

ReDim count(ng) As Integer

For i=1 To a

mid=Int(Z2(i)/d)

an=360/(2*(mid+1)1)

mid=mid^2+Int(Z1(i)/an)+1

count(mid)=1

Next i

For i=1 To ng

If (count(i)=1) Then

r=r+1

End If

Next i

Cells(4,3)=r

End Sub

3实例应用

大岗山水电站是大渡河干流的大型水电工程之一。坝区节理裂隙发育广泛，主要发育7组优势裂隙。本文以坝区左岸平硐、建基面①至⑥组裂隙资料为例进行实证分析。图5给出了左岸坝区建基岩体开挖前后第①组节理产状极点分布的网格剖分图，其它各组与其相似。

依据上述节理产状极点分布分形维D的计算方法和编程来计算左岸坝区建基岩体开挖前后各组裂隙产状极点分布的分形维D。整个程序仅涉及一个输入变量，即环数n。依据数据统计可知，当n=1，2，3等较小值时，极点占有的小方格数N(δ)较少，很难满足分形要求；而当n>45等较大值时，极点占有的小方格数N(δ)基本相同，同样不具分形特性。所以，将n值设在6~45之间，取整数，相应的无标度区间约为（ln4，ln8）（图6(a)6(b)）。

图5n=15时左岸第①组裂隙产状极点图网格剖分

Fig.5Schmidt subdivision of ① joint orientation

at left dam foundation when n=15

图6节理产状极点分布对应的lnδ~lnNδ曲线

Fig.6Curves of lnδ~lnNδ corresponding to the

orientation pole distribution for joints

表3为左岸坝区平硐、建基面各组节理产状极点分布的分形维D计算结果。可以看出，分形维D可以较好地阐释建基面开挖前后节理产状分布的差异性，其大小不仅与节理数有关，更重要的是与节理产状分布的分散度有很大的关系。具体分析如下。

表3左岸坝区平硐、建基面各组节理产状极点分形维

Table 3Fractal dimension of orientation pole distribution

for joints collecting from adit and foundation surface in left dam

节理组平硐建基面节理数/条分形维D相关系数节理数/条分形维D相关系数①2810.5060.9926090.5370.991②5470.4790.9945130.5960.992③3400.2690.9945680.4920.971④7240.5270.99622230.5750.994⑤1450.2080.9892160.3550.978⑥1 0200.7260.99411670.5930.992①至⑥3 0570.8070.99252960.8140.985（1）一般情况下，节理数越多，产状极点分布的分形维D越大。同组节理，依据建基面测网资料计算的节理产状极点分布分形维D较平硐中的大，主要由于随着建基面的开挖节理条数相对增多，节理产状信息更加丰富，如第①、③、④、⑤组节理，随着节理条数增大，分形维D明显增大。非同组节理间也遵循该规律，如第④、⑥组节理条数大于第①、②、③、⑤组节理，其分形维D也相应较大。

（2）分形维D可用来描述节理产状极点分布的分散程度，分散度越大，分形维D越大。一般情况，节理数越多分形维D越大，同时意味着节理产状极点分布地较为分散；如果节理数较小，而分形维D较大，则说明在节理产状极点图中相同倾向倾角的产状较少，极点重合率较小，因此增大了产状极点分布的范围。如表3中的第①、②组节理，建基面开挖后，第②组节理的条数小于第①组节理，但是其分形维D却相对较大，说明第②组节理产状的分散度较大。

（3）节理产状极点分布的分形维D在理论上满足条件：0≤D＜1［7,16］。如果某组节理中各条节理的产状完全相同，则所有节理的极点在Schmidt等面积投影网中将重合在1个极点上，该情况下，无论怎样改变投影网中小方格的边长δ，被极点占有的小方格数N(δ)恒等于1，而lnδ~lnN(δ)曲线必是一条水平直线，其斜率恒等于0，即说明在这种理想状况下，节理产状分形维D=0。在欧氏几何空间中，1个点的维数为0，所以该理想状况也可以用欧氏几何中对点的维数来解释。实际上，理想状态一般不存在，即使同组裂隙中各条节理产状也不尽相同，因此其极点分布具有一定的范围，所以一般情况下，节理产状分形维D＞0。如果某一数量的节理，其中各条节理产状差异较大，则极点图中极点的分布范围较广，即产状的分散度较大，则分形维D较大，这与计算结果中的第（2）条结论相一致。另外，由未分组的①至⑥组节理产状极点的分形维可知，节理产状极点分布分形维的值恒小于1。因为只有当极点占有小方格数N(δ)与每次循环剖分的小方格总数N相等时，lnδ~lnN(δ)曲线的斜率才等于1。然而，由于点本身的几何特性，点与点之间总是存在一定的空间（或者说，极点是无法连续地充满Schmidt图的），所以当无限度地对Schmidt图进行网格剖分时，N(r)＜N恒久成立。所以，即使节理数很多，产状分散度很大，各个方向的节理均有分布，其产状分形维D也恒小于1。

4结论

（1）本文针对节理产状极点分布分形维D的计算原理，进行了程序设计。在设计过程中，利用极坐标将节理产状极点进行数值化处理，同时通过将小方格极坐标域化，可直接以被相应数值化的极点占有的极坐标域的个数来代替被极点占有的小方格个数，使编程更为简便。

（2）采用Visual Basic for Applications（VBA）进行编程，在Excel中直接进行数据处理，使程序界面操作简单，容易实现。

（3）根据对大岗山水电站坝区左岸6组优势节理产状极点分布的分形维D的计算，验证了程序的正确性与可用性。

参考文献（References）：

［1］石根华.岩体稳定分析的赤平投影方法［J］.中国科学,1977(3):260271.(SHI Genhua,Stereographic Projection method of Rock mass Stabilization in Rock mass Stability Analysis［J］.Science in China,1977(3):260271.(in Chinese))

［2］何绍勋.赤平投影在构造地质学中应用的进展［J］.地质科技情报,1989,8(3):3543.(HE Shaoxun.Application of Stereographic Projection in Structural Geology［J］.Geological Science and Technology Information,1989,8(3):3543.(in Chinese))

［3］吴绍强.极射赤平投影法在岩质边坡稳定性分析中的应用［J］.西部探矿工程,2009,(10):117121.(WU Shaoqiang.Application of Stereographic Projection Method in Rock Slope Stability Analysis［J］.WestChina Exploration Engineering,2009,(10):117121.(in Chinese))

［4］高玉娟,徐付玲,张亚明.赤平投影方法在大地构造研究中的改进及应用［J］.辽宁工程技术大学学报(自然科学版),2009,28(S1):6567.(GAO Yujuan,XU Fuling,ZHANG Yaming.Improvement and Application of Stereographic Projection in Tectonic［J］.Journal of Liaoning Technical University(Natural Science),2009,28(S1):6567.(in Chinese))

［5］王新刚,陈晋,胡斌,等.基于赤平投影的节理控制性黄土边坡破坏模式分析［J］.煤炭安全,2013,44(8):210212.(WANG Xingang,CHEN Jin,HU Bin,et al.Failure Mode Analysis of the Jointed Controlled Loess Slope Based on the Stereographic Projection［J］.Safety in Coal Mines,2013,44(8):210212.(in Chinese))

［6］刘梦君,许模,隆振锐.某水电站厂房调压室半确定性块体稳定性浅析［J］.南水北调与水利科技,2011,9(6):136139.(LIU Mengjun,XU Mo,LONG Zhenrui.Stability Analysis of SemiDeterministic Blocks for the Surge Chamber in the Factory Building of a Hydropower Station［J］.Southtonorth Water Diversion and Water Science & Technology,2011,9(6):136139.(in Chinese))

［7］陈剑平,王清,谷宪民,等.岩体节理产状极点分布的分形维［J］.岩石力学与工程学报,2007,26(3):501508.(CHEN Jianping,WANG Qing,GU Xianmin,et al.Fractal Dimension of Orientation Pole Distribution for Rock Mass Joints［J］.Chinese Journal of Rock Mechanics and Engineering,2007,26(3):501508.(in Chinese))

［8］FISHER R.Dispersion on a sphere［J］.Proceedings of the Royal Society of London,Series A,1953,217:295305.

［9］BINGHAM C.Distribution on the sphere and on the Projective Plane［Ph.D.Thesis］［D］.New Haven:Yale University,1964.

［10］GOODMAN R E.不连续岩体中的地质工程方法［M］.北京:中国铁道出版社,1980:3381.(GOODMAN R E.Methods of Geological Engineering in Discontinuous Rock Mass［M］.Beijing:China Railway Press,1980:3381.(in Chinese))

［11］陶振宇,潘别桐.岩石力学原理与方法［M］.武汉:中国地质大学出版社,1981:2330.(Tao Zhenyu,Pan Bietong.Rock Mechanics Principles and Methods［M］.Wuhai:China University of Geosciences press,1981:2330.(in Chinese))

［12］范雷,王亮清,唐辉明.节理岩体结构面产状的动态聚类分析［J］.岩土力学,2007,28(11):24052408.(FAN Lei,WANG Liangqing,TANG Huiming.Dynamic Cluster Analysis of Discontinuity Orientations of Jointed Rock Mass［J］.Rock and Soil Mechanics,2007,28(11):24052408.(in Chinese))

［13］孙宪春,万力,蒋小伟.节理产状分组的k均值聚类分析及其分组结果的费歇尔分布验证法［J］.岩土力学,2008,29(S):533537.(SUN Xianchun,WAN Li,JIANG Xiaowei.Effective Categorization of Joints by kmeans and Cluster Analysis Its Verification by Fisher Distribution ［J］.Rock and Soil Mechanics,2008,29(S):533537.(in Chinese))

［14］冯羽,马凤山,巩城城,等.节理岩体结构面优势产状确定方法研究［J］.工程地质学报,2011,19(6):887892.(FENG Yu,MA Fengshan,GONG Chengcheng.et al.Data Analysis Method for Optimized and Dominant Orientations of Joints in Rock Mass［J］.Journal of Engineering Geology,2011,19(6):887892.(in Chinese))

［15］隆振锐,许模,张翔,等.基于概率模型的随机裂隙优势产状取值［J］.南水北调与水利科技,2012,10(1):121124.(LONG Zhenrui,XU Mo,ZHANG Xiang,et al.Discussion of Values Dominant Orientations of Random Fractures Based on Probability Models［J］.Southtonorth Water Diversion and Water Science & Technology,2012,10(1):121124.(in Chinese))

［16］宋丽娟,许模,卢书强,等.节理产状极点分布的分形维改进方法［J］.岩石力学与工程学报,2013,32(S2):33033308.(SONG Lijuan,XU Mo,LU Shuqiang,et al.Fractal Dimension Improved Algorithm of Orientation Pole Distribution for Joints［J］.Chinese Journal of Rock Mechanics and Engineering,2013,32(S2):33033308.(in Chinese))

推荐访问:产状 节理 极点 计算方法 程序