基于体视学技术的二维数据模拟三维裂隙网络

倪春中，刘春学，张世涛，李爱敏，李雨健，苍宝海，刘 雷

(1.昆明理工大学 国土资源工程学院，云南 昆明 650093; 2.云南财经大学 城市与环境学院，云南 昆明 650221;3.云南锡业集团有限责任公司，云南 个旧 661021; 4.中国石油 辽河油田分公司 特种油开发公司，辽宁 盘锦 124010;5.中国地震局 地震预测研究所，北京 100038)

基于体视学技术的二维数据模拟三维裂隙网络

倪春中1，刘春学2，张世涛1，李爱敏3，李雨健1，苍宝海4，刘 雷5

(1.昆明理工大学 国土资源工程学院，云南 昆明 650093; 2.云南财经大学 城市与环境学院，云南 昆明 650221;3.云南锡业集团有限责任公司，云南 个旧 661021; 4.中国石油 辽河油田分公司 特种油开发公司，辽宁 盘锦 124010;5.中国地震局 地震预测研究所，北京 100038)

CT扫描；体视学；裂隙网络；岩石裂隙；个旧锡矿

近年来国内外对于裂隙空间分布的研究虽然很多，可所用的数据大多是通过钻孔或者露头等获取的一维或二维样品资料，而真正的三维裂隙网络分布样品资料是很难观测的，因此不少学者尝试通过统计学和应用概率等方法将一维或二维的数据推算三维岩石裂隙网络的属性。地下岩石裂隙的勘察方法主要有钻孔岩心法、超声成像法、钻孔电视法以及钻孔流量测井法等[1]，小尺度的观测方有如核磁共振(NMR)、X-射线扫描(CT)、γ-射线扫描和同步加速器(synchrotron)等。裂隙空间分布具有很多属性，各属性之间大部分不是相互独立的，存在多种形式的相关关系，因此确定裂隙的长度、密度等重要的属性后[2-3]，其他变量可以随之确定。Berkowitz和Adler发现裂隙网络与观测截面交嵌的圆盘直径以及其他变量均服从Power Law分布、lognormal Law分布和Exponential Law分布[4]。Adler和Thovert运用二维图像中的痕迹长度及数量估算三维裂隙网络中裂隙密度和概率密度等属性[5]。Yeomans认为裂隙的个别属性(如长度、位移、开口宽度等)以及整个裂隙网络分布也服从Gamma Law趋势[6]，而Dershowitz和Einstein则认为裂隙方向可以用Bivariate Normal分布和Fisher分布进行描述[7]。

本研究以云南个旧锡矿区高松矿田为例，所用的小尺度裂隙数据取自样本岩石的CT扫描图像，应用体视学技术的方法将二维图像中的裂隙延展及其连接等数据反演到三维空间，然后在三维空间相关特征的约束条件下，反演推断出三维空间中裂隙密度分布，最终实现应用二维裂隙数据模拟三维裂隙网络。

1 地质概况

高松矿田是云锡公司的主要生产矿山之一，位于个旧锡矿东矿区北部，处于矿区的五子山复背斜北段，

夹持于南北向个旧断裂、甲界山断裂与东西向个松断裂、背阴山断裂之间，地质条件复杂。高松矿区主要出露地层为二叠系和三叠系，以白云岩、白云质灰岩为主，另有古近系与新近系泥岩、第四系残坡积物零星出露。

褶皱构造有近EW向对门山-阿西寨向斜、近SN向五子山背斜等。断裂构造可以分为EW，NE，NW，SN四组，其中以近EW向和NE向为主，NW向次之。EW向断裂组自北而南有个松断裂、麒麟山断裂、马吃水断裂、高阿断裂、背阴山断裂，呈近等距离分布，其中个松和背阴山断裂规模较大，麒麟山断裂次之。NE向断裂组自西而东有莲花山断裂、芦塘坝断裂、麒阿西断裂。NW向断裂组主要为大箐东断裂、黑码石断裂、驼峰山断裂、阿西寨断裂。SN向断裂组不发育(图1)。

2 样本裂隙提取

2.1 岩石样本

本研究在个旧锡矿高松矿田1 360，1 540，1 720和1 920 m四个不同中段采取岩石样品35件，从中选取了12件岩石样品进行CT扫描，用于微观尺度下三维裂隙网络的模拟，各样品采样平面位置如图1所示，在中段的位置分布如表1。

图1 个旧矿区高松矿田构造地质略图(文献[8]，修改)Fig.1 Tectonic geological sketch of Gaosong mine in GejiuT2g1.中三叠统个旧组卡房段白云岩灰岩；T2g2.中三叠统个旧组马拉格段白云岩；T2g3.中三叠统个旧组白泥硐段灰岩

样品编号001002003004005006007008009010011012采样位置/m192019201720172015401540154015401360136013601360

2.2 裂隙迹线提取

CT图像中提取的信息主要包括：扫描定位的上下、前后、左右以及图像中特殊点的位置，扫描定位后，CT机给出设定的图像放大倍数，图像中以X/Y的数值给出光标指示点的二维坐标；试样的几何特征，包括试验中试样直径、长度、两点间的距离和角度[9-10]，进行几何测量和内部变形的研究；图像的分布特征，从CT数及其方差正负变化的组合情况，可以判断试件整体或某个部位处于压密还是裂隙扩展等；岩土试样中某些节理构造，在CT图像中表现为某些方向的明暗条纹。

本次试验按照岩石的实际产状对岩样进行CT扫描，得到了岩样的三维空间网格点上的CT数，用Vo-xelCalc NDT Software进行了初步数据处理，得到[xyzCTN]格式的元数据。然后用Matlab编写了数据处理程序，对CT扫描得到的元数据进行了进一步处理。通过调整各个岩样CT数的阈值，使CT数在空间的变化更加明晰，便于观察裂隙的展布(图2)。

从50个水平扫描剖面中提取了相关的裂隙数据。裂隙数据除了实际观测的迹线外，另外一部分资料利用Matlab编写的辅助程序从CT扫描图像。岩样水平剖面的大小多数为60 mm×80 mm。从中提取了284条裂隙迹线，其中最长为64.1 mm，最短为2.2 mm，平均值为16.69 mm，方差为12.33 mm，集中在10 mm左右。

3 基于体视学方法由二维迹线模拟裂隙长度及密度

体视学利用形态学和数学的分析，借助计算机及数据处理系统和显微镜及显做成像系统，将二维平面

经过成像及计算机分析处理得到三维形态，准确地对物体进行定量及形态结构分析。本文应用体视学技术，研究圆盘裂隙三维空间与二维剖面之间的相互联系，解析裂隙长度在不同维数之间的联系和转化，以处理二维剖面中的样本裂隙迹线数据。

定义直径为φ的圆盘，其中心点与观测面之间的距离为z，其单位法线矢量与平面法线之间的夹角为α(图3)。

假设φ为常数，直径为φ的圆盘与观测面交切的比例为S(z,φ)。从图3可以看出，只有当|z|<φ/2时才会有交切。对于给定的z，S(z,φ)可以表示为

(1)

式中：z为裂隙圆盘中心到观测面S的距离，mm；φ为圆盘直径，mm。

下述公式应用幂次定律可以计算裂隙直径分布[11-13]：

(2)

式中：α为分形维数，对于裂隙迹线长度来说，1≤α≤3；a为待定参数，由式(3)确定。

如果裂隙直径长度假定在φm≤φ≤φM的范围内，那么常数α可以表示为：

(3)

(4)

式中：φM,φm分别为裂隙最大直径和最小直径，mm；α为分形维数，根据裂隙分布试验确定。

根据直径的幂次定律，假定裂隙直径最大值φM是无限的，迹线c的概率可表示为：

图2 岩样CT图像裂纹提取Fig.2 Crack extraction from rock CT image

图3 利用体视学技术研究圆盘裂隙三维空间与二维剖面之间相互关系原理Fig.3 Study on the relationship between 3D space and 2D profiles of disc fissures used Stereologya.直径为φ的裂纹圆盘中心点与观测面S的距离为z；b.裂纹圆盘与观测面S的交切示意图

(5)

(6)

(7)

由于迹线长度总小于或等于裂隙直径，为了便于计算，设定c=φ，则有：

(8)

因为c1≤φ≤φM,c0=φM时h(φ)=h1，可得：

此次从50个水平面中提取的裂隙迹线长度的频率直方图见图4。可见岩样剖面中的裂隙迹线长度分布遵从对数正态分布，与一般研究中的发现一致。在考虑取样偏差后，为了研究的简便，可以认为迹线长度分布服从幂律分布，即将较短的裂隙迹线长度省去，与一般的裂隙观察一致，从而可以用前述的理论方法进行分析。

在运用公式(9)计算裂隙直径时，依据实际观测数据，最小直径φm设定为1 mm；最大直径φM在1 mm到200 mm之间变化，本文取最大值200 mm；指数α一般选择两个值，分别是1.5和2.5，根据与现场裂隙分布相比较，取α=2.5，可以发现该值与实际情况相吻合。依据公式(2)—公式(6)，可分别求出α=1.500 5，β=0.422 6，b=1.5。

图4 岩样裂隙迹线长度分布频率Fig.4 Frequency distribution of trace length of rock fissures

4 结论

表2 裂纹直径和密度计算结果

致谢：在野外工作期间，得到云南锡业公司松矿朱文捷高工的大力帮助，在此表示感谢！

[1] 张连明，周峰，张兆坤.地下岩层裂隙的几种勘察研究方法[J].水文地质工程地质，2003，4：111-114. Zhang Lianming,Zhou Feng,Zhang Zhaokun.Several methods for investigating crack in underground rock stratum[J].Hydrogeology and Engineering Geology，2003,4:111-114.

[2] 白斌，邹才能，朱如凯,等.四川盆地九龙山构造须二段致密砂岩储层裂缝特征、形成时期与主控因素[J].石油与天然气地质，2012，33(4)：526-535. Bai Bin，Zou Caineng，Zhu Rukai,et al.Characteristics,timing and controlling factors of structural fractures in tight sandstones of the 2nd member of Xujiahe Formation in Jiulong Mountain structure,Sichuan Basin[J].Oil & Gas Geology，2012，33(4)：526-535.

[3] 张仲培,刘士林,杨子玉,等.塔里木盆地麦盖提斜坡构造演化及油气地质意义[J].石油与天然气地质,2011,32(6):909-919. Zhang Zhongpei,Liu Shilin,Yang Ziyu,et al.Tectonic evolution and its petroleum geological significances of the Maigaiti Slope,Tarim Basin[J].Oil & Gas Geology，2011,32(5):909-919.

[4] Berkowitz B,Adler P M.Stereological analysis of fracture network structure in geological formations[J].Geophysics Research,1998，B103:15339-15360.

[5] Adler P M，Thovert.Fractures and fracture network[M].USA:Klu-wer Academic Publishers,1999.

[6] Acuna J,Yeomans Y.Application of fractal geometry to the study of networks of fractures and their pressure transient[J].Water Resources

Research,1995,31(3):527-540.

[7] Dershowitz W S，Einstein H H.Characterizing rock joint geometry with joint system models[J].Rock Mechanics Rock Engineering,1988,21:21-51.

[8] 庄永秋，王任重，郑树培，等.云南个旧锡铜多金属矿床[M].北京：地震出版社，1996. Zhuang Yongqiu,Wang Renzhong,Zheng Shupei,et al.Geology of Gejiu tin-copper polymetallic depositl(in China)[M].Beijing:Seismological Press,1996.

[9] 苏建政，黄志文，龙秋莲，等.基于ANSYS软件的降低破裂压力机理模拟[J].石油与天然气地质，2012，33(4)：640-645. Su Jianzheng,Huang Zhiwen,Long Qiulian,et al.ANSYS-based si-mulation of fracturing pressure reducing mechanism[J].Oil & Gas Geology,2012，33(4)：640-645.

[10] 苏娜，段永刚，于春生.微CT扫描重建低渗气藏微观孔隙结构——以新场气田上沙溪庙组储层为例[J].石油与天然气地质,2011，32(5)：792-796. Su Na,Duan Yonggang,Yu Chunsheng.Reconstruction of microsco-pic pore structure in low permeability gas reservoir by micro-CT scanning:an example from the Upper Shaximiao Formation in Xinchang gas field[J].Oil & Gas Geology，2011，32(5)：792-796.

[11] Barton C A，Zoback M D.Self-similar distribution and properties of macroscopic fractures at depth in crystalline rocks in the Cajon Pass scientific drill hole[J].International Joural of Rock Mechanics and Ming Science，1993，30(1)：5181-5200.

[13] 倪春中，刘春学，张世涛.从岩石露头裂隙迹线估算裂隙三维空间方向[J].石油与天然气地质，2013，34(1)：102-106. Ni Chunzhong,Liu Chunxue,Zhang Shitao.Estimation of three-dimensional distribution of fissures according to fissure traces on outcrops[J].Oil & Gas Geology,2013,34(1)：102-106.

(编辑 董 立)

Stereology-basedsimulationof3Dfracturenetworkwith2Ddata

Ni Chunzhong1，Liu Chunxue2，Zhang Shitao1，Li Aimin3，Li Yujian1，Cang Baohai4，Liu Lei5

(1.EngineeringInstituteofLandandResources,KunmingUniversityofScienceandTechnology,Kunming,Yunnan650093,China；2.SchoolofUrbanandEnvironment,YunnanUniversityofFinanceandEconomics,Kunming,Yunnan650221,China；3.YunnanTinGroupCo.Ltd.,GejiuYunnan661021,China；4.SpecialOilDevelopmentCompany,PetroChinaLiaoheOilfieldCompany,Panjin,Liaoning124010,China; 5.InstituteofEarthquakePrediction,ChinaSeismologicalBureau,Beijing100038,China)

Various attribute information of fissures such as length,width and density can be extracted from CT scanning images of oriented rock samples that are systematically collected through field work.These information can play a important role in defining spatial distribution of fissure networks within rocks.Both trace distribution and fissure network distribution follow the power-law distribution,thus statistical analysis can be performed on the field-collected trace data.The extension length and interlink of fissures on 2D images can be inversed to 3D space by applying stereological techniques,so as to deduce 3D fissure length and density distribution and realize conversion of rock fissure attributes from 2D to 3D.This approach was applied in Gaosong ore field.Based on analysis of the traces on the CT scanning images,we obtained the maximum trace length of 64.1 mm,the minimum of 2.2 mm,a mean of 16.69 mm and a variance of 12.33 mm.The fissures were divided into 9 groups,and the probability distribution of diameters of each group was obtained.Simulation was performed on the extracted trace data when the fractal dimension(α)is 1.5 and 2.5 respectively.The mean of fissure diameter was 27.60 mm,and fracture network density was 0.30/cm3.This approach can be popularized to other mines to obtain attributes of 3D fissure network when only 2D data are available.

CT section image,stereology technology,fissure network,rock fissure,Gejiu tin deposit

2013-07-02;

：2013-12-01。

倪春中(1979—)，男，博士、讲师，构造地质学。

国家自然科学基金项目(40902058);云南省应用基础研究计划项目(2013FZ026);云南省教育厅科学研究基金项目(2011Y363)。

0253-9985(2014)01-0148-05

10.11743/ogg20140119

TE19

：A