从岩石露头裂隙迹线估算裂隙三维空间方向

倪春中，刘春学，张世涛

(1.昆明理工大学国土资源工程学院，云南昆明650093; 2. 云南财经大学城市与环境学院，云南昆明650221)

裂隙是地球科学中的一类重要构造，在地壳中广泛分布，与裂隙网络渗流关系密切，因此查清裂隙的空间分布显得尤为重要。但目前获取裂隙空间分布特征的途径，除了在实验室内进行小规模的分形研究外［1-4］，大多数是通过钻孔，坑道开挖面以及一些岩石露头获取。野外裂隙调查的方法主要有测线法［5-7］和统计窗法［8］。但在实际工作中存在两个问题:首先是不同的取样方法，人为测量偏差以及仪器误差都有可能影响结果精度;其次在野外由于岩石松动以及岩块掉落等原因，运用上述两种方法直接对陡峭的岩墙进行测量存在很大的安全隐患。本文通过对数字影像提取裂隙痕迹，克服了上述方法的缺点。利用数字影像处理软件，提取裂隙痕迹属性，包括走向、倾角、裂隙间距、大小以及形状等。裂隙三维空间分布与裂隙痕迹的这些属性有关，尤其是走向、倾角。每组裂隙的走向、倾角都有一个均值，围绕着这个均值散布着众多裂隙。一般认为裂隙面法线方向的概率密度可以用Fisher 分布进行描述，运用Fisher分布密度函数的前提条件是对裂隙进行分组，根据参数K，模拟出裂隙面产状，同时由于该函数是可积函数，产生随机数较方便。本文先考虑已知岩墙面的产状，裂隙平均走向、倾角以及Fisher 分布参数K，计算岩墙面裂隙痕迹与水平面的夹角，估算出裂隙网络面的属性参数，在此基础上，利用模拟退火法对原问题求逆，并将该方法应用于个旧锡矿区高松矿田，考察其预测裂隙网络分布的有效性。

1 裂隙面产状确定其迹线的射线角

如前所述，裂隙方向对裂隙的其它属性影响至关重要。每组裂隙都有一个均值，裂隙方向围绕着这个均值分布，Fisher 分布参数K 可以用来生成这些离散值。

1.1 Fisher 分布

岩体裂隙面在空间的倾角用其法线的方向表示。应用球面坐标系可以给出Fisher 分布模型的表达式［9-10］为:

式中:θ＇为经过变换后的裂隙面倾角，(°);K 表示裂隙组极点集的集中程度参数，无量纲。对岩石裂隙来说，K 值取值一般介于20～300，其取值方法参见文献［5］。根据该公式可以推导出Fisher分布函数在极点图上用极坐标表示的面积密度，即以真实平均值θ 为中心的角面积dA 中出现裂隙的概率，其函数表示为:

式中:θ 为裂隙面倾角真实平均值，(°)。当θ =0时，f(θ)值最大。该式反映围绕裂隙均值的对称分布。K 值越大，方向变量越集中于真实平均值。

可以用下述公式近似估计计算机产生的Fisher分布［11］:

式中:Δθ 为裂隙与均值的差，(°);Random(0，1)为介于0 到1 的均匀随机数。为了确定新的倾角和倾向，Δθ 围绕平均值旋转一个介于0～360°的任意角。

1.2 裂隙方向和迹线角之间的关系

定义节理平面法向量的倾向和倾角为(αj，βj)，岩石面法向量的倾向和倾角为(αf，βf)(图1)。

首先根据岩石面的产状(αj，βj)可以分别计算垂直于节理平面(j 向量)和岩石面(f 向量)的法线坐标，迹线矢量t 为矢量j 与矢量f 的叉乘t =j×f，裂隙走向的单位矢量可以表示为s =(sx，sy，sz)=(cosαf，sinαf，0)。

裂隙迹线角θt可以用s 与t 的点乘确定:

式中:θt为迹线与水平面的锐夹角，(°);sx，sy，sz分别为裂隙走向单位矢量在x，y，z 轴上的分量;tx，ty，tz分别为迹线方向单位矢量在x，y，z 轴上的分量。

图1 裂隙方向与迹线角关系Fig.1 Relationship between fissure direction and trace dip

由于θt可出现在不同的象限内，不是岩石面上迹线角的唯一值，为了确保迹线角的唯一性，本文采用射线角θrake表示在岩墙面中从水平面到迹线的顺时针角，它可作为迹线角的唯一度量值。

2 二维裂隙痕迹反演其三维空间分布

基于以上数据，本文利用模拟退火法，对逆问题求解，就可以得到裂隙的倾向和倾角。

2.1 模拟退火法

模拟退火优化方法最早由Kirkpatrick 和Cerny 等人提出和研究［12-13］，分析了优化领域里的许多优化组合和运用统计物理求解物理系统最低能量状态问题，并基于此导出了模拟退火寻优法的基本框架。

考虑一个物理系统s，设s 共有N 个状态。当达到平衡后，系统将随机地处于一种状态，并且处于状态xi的概率满足Boltzmann 分布，即:

式中:ui为s 处于状态xi时具有的能量，J;Tk为系统处于热平衡状态下对应的温度，℃。

从Boltzmann 分布可以看出，温度越低，系统处于具有小能量状态的概率也越大，当温度趋于绝对零度时系统处于具有最小能量状态的概率趋于1。所以，为了使物体s 具有理想的规则晶格结构(对应于能量最小的状态)，通常的作法是给s加足够高的温度，使其达到液化状态，然后缓慢地逐步降低温度，直到温度接近于绝对温度。这样，s 将随着温度的下降，逐渐地从高能量状态过渡到能量越来越低的状态，直至能量小或接近最小。这个过程，物理上称为退火过程。

显然，如果把优化问题的状态看成物理系统的状态，目标函数值看成物理系统的能量，则优化问题和物理系统求得最低或接近最低能量状态问题是一致的。模拟退火寻优法就是基于此而导出的。

2.2 求解裂隙面产状

由统计物理的退火过程，可以构造如下模拟退火仿真优化算法。

1)设定初始状态X(0)，定义节理平面法向量的倾向为αj和倾角为βj，并假定αj0=360°和倾角βj0=90°，C=C0，C0为设置的控制参数。

2)随机产生新状态X(m)，倾向αjm和倾角为βjm，估计优化目标值:

式中:μ，σ，θ 分别表示平均值、标准偏差以及偏斜度;A，B，C 为常数;I 表示二维图像中痕迹方向角分布的参数;T 表示根据公式(2)-公式(9)由模拟退火法产生的参数。

3)计算转移概论P:

式中:Xm+1为迭代计算产生的新解;f(Xm+1)为式(6)代入新解求得的优化目标函数值;c 为初始控制参数。

4)判断接受或放弃新解Xm+1。当P(Xm+1=Xm)＞r(r 为设定的允许误差)，接受Xm为新解;否则，放弃Xm+1为新解。

5)判断终止条件，满足条件则停止迭代;否则转向2)。

3 实例分析

为了考察二维裂隙数据估算裂隙三维空间分布的有效性［14-15］，利用个旧锡矿区高松矿田某一开挖面上得到的裂隙迹线数据进行了实例分析，将计算得到的节理产状与实测数据进行了比较。图2a 为个旧锡矿区高松矿田某临矿山公路的开挖面，长16.2 m，高约3 m。该岩墙面的产状为255°∠80°，岩性主要为白云质灰岩，发育有多组节理。经数字影像处理，获取179 条裂隙痕迹，提取了裂隙痕迹倾角三维分布图(图2b)。根据裂隙痕迹与水平面夹角θrake值分布直方图，岩墙面裂隙痕迹以90°为界分为两组(图3)，均值分别为49.0°和128.2°。

图2 岩墙面裂隙痕迹分布影像a)相对应的裂隙痕迹三维分布b)Fig.2 Distribution of fissure trace in rock surface image (a)corresponding with the 3D distribution of the fissure trace (b)

图3 裂隙痕迹θrake值频度直方图Fig.3 Histogram showing the value of fissure trace

依据裂隙丛聚规律，节理面倾向取值范围为［0°，359°］，倾角取值范围为［0°，90°］。根据公式(2)，计算出K 值分别等于20，50，100，150 和200 时Fisher 分布。图4 表示符合Fisher 分布，倾向和倾角均值分别为120°和20°的节理组，图4a 中k=50，图4b 中k=200。本文通过与现场比较，取K 值等于200，计算中所用的参数A，B，C 所占的权重分别取70%，30%和0，求得误差5%范围内θrake和SD(标准方差)值(表1)。

为了求得裂隙面的三维空间分布，根据所提取裂隙痕迹数据，利用Matlab 语言编程，然后根据最优函数公式(6)及判别式(7)求解最小误差(设定r值≤0.05)。裂隙面产状估算结果如表1 所示，可见模拟情况比较接近岩体裂隙的实际分布情况。

4 结论

图4 符合Fisher 分布的节理组极点图(倾向120°，倾角20°)Fig.4 Point diagram of joint set in agreement with the Fisher distribution (mean dip direction of 120°and mean dip of 20°)

表1 裂隙痕迹模拟结果Table 1 Simulation of fissure traces

由岩墙面裂隙痕迹估算裂隙的空间产状，可以通过影像数字处理，运用计算机自动提取裂隙迹线信息，有效克服了在陡峭的岩墙面上无法运用测线法和统计窗法进行裂隙测量的缺点，消除了人为运用上述两种方法造成的误差，同时也可以弥补在野外工作中遗漏的部分裂隙信息。裂隙迹线进行分组后，确定每组裂隙迹线倾角均值，根据现场裂隙分布的集中程度，可以确定Fisher 分布常数K，结合岩墙面方向以及镜头方向与岩墙面法线夹角，利用模拟退火优化仿真算法，估算得到裂隙的倾向和倾角，实现从二维裂隙生成三维裂隙网络，大大减少了人工现场测试裂隙工作量。

［1］陈田勇，毛鑫，刘仕银，等. 利用分形理论计算相对渗透率曲线——以南襄盆地双河油田核桃园组六油组为例［J］.石油与天然气地质，2012，33(4):578-581.Chen Tianyong，Mao Xin，Liu Shiyin. Fractal theory-based calculation method of relative permeability curves—a case study from the Hetaoyuan Formation in Shuanghe oilfield，Nanxiang Basin［J］.Oil ＆ Gas Geology.2012，33(4):578-581.

［2］王津义，彭金宁，王彦清，等.贵州南部印支晚期——燕山期构造变形特征［J］.石油实验地质，2012，34(4):362-367.Wang Jinyi，Peng Jinning，Wang Yanqing，et al.Tectonic deformation features from Late Indosinian to Yanshanian in southern Guizhou［J］.Petroleum Geology ＆ Experiment，2012，34(4):362-367.

［3］冯阵东，戴俊生，邓航，等. 利用分形几何定量评价克拉2气田裂缝［J］.石油与天然气地质，2011，32(6):928-933.Feng Zhendong，Dai Junsheng，Deng Hang，et al. Quantitative evaluation of fractures with fractal geometry in Kela-2 gas field［J］.Oil ＆ Gas Geology.2011，32(5):928-933.

［4］郭恺，仝兆岐，滕厚华.地震采集面元尺度与成像横向分辨能力分析［J］.石油与天然气地质，2012，33(1):141-147.Guo Kai，Tong Zhaoqi，Teng Houhua. Analysis on bin size of seismic survey and lateral resolution of imaging［J］.Oil ＆ Gas Geology.2012，33(1):141-147.

［5］Priest S D，Hudson J A.Estimation of discontinuity spacing and trace length using scanline surveys［J］.International Journal of Rock Mechanics，Mining Sciences and Geomechanics Abstracts，1981，18(3):183-197.

［6］Cruder D M.Describing the size of discontinuities［J］.International Journal of Rock Mechanics and Mining Sciences and GeomechanicsAbstracts，1977，l4(3):133-137.

［7］余一欣，周心怀，徐长贵，等. 渤海海域新生代断裂发育特征及形成机制［J］.石油与天然气地质，2011，32(2):273-279.Yu Yixin，Zhou Xinhuai，Xu Changgui，et al. Characteristics and formation mechanism of the Cenozoic faults in the Bohai Sea waters［J］.Oil ＆ Gas Geology，2009，30(2):273-279.

［8］Paul P H. Estimating the mean length of discontinuity trace［J］.International Journal of Rock Mechanics and Mining Sciences and GeomechanicsAbstracts，1981，18(3):221-228.

［9］陆峰，王俊奇.Fisher 模型在岩体裂隙面模拟中的应用［J］.中国水利水电科学研究院学报，2011，8(4):309-313.Lu Feng，Wang Junqi.Application of Fisher Model in rock fracture simulation［J］. Journal of China Institute of Water Resources and Hydropower Research，2011，8(4):309-313.

［10］Liu Chunxue，Koike K. Extending multivariate space-time geostatistics for environmental data analysis［J］.Mathematical Geology，2007，39(2):289-305.

［11］Priest S D. Discontinuity analysis for rock engineering［M］.London:Chapman ＆ Hall，1993:473.

［12］Kirkpatrick S，Gelatt C D，Vecchi M P.Optimization by simulated annealing［J］.Science，1983，220(4598):671-680.

［13］Cerny V.Thermodynamical approach to the traveling salesman problem:an efficient simulation algorithm［J］.Journal of Optimization Theory and Applications，1985，45(1):41-51.

［14］Koike K，Liu Chunxue，Sanga T. Incorporation of fracture directions into 3D geostatistical methods for a rock fracture system［J］. Environment Earth Science，2012，66(5):1403-1414.

［15］李旭兵，刘安，曾雄伟，等.雪峰山西侧地区寒武系娄山关组碳酸盐岩储层特征研究［J］. 石油实验地质，2012，34(2):153-157.Li Xubing，Liu An，Zeng Xiongwei，et al. Characteristics of carbonate reservoirs of Cambrian Loushanguan Formation to the west of Xuefeng Mountain［J］. Petroleum Geology ＆ Experiment，2012，34(2):153-157.