利用地质统计学方法模拟岩石裂隙网络的三维空间分布——以云南个旧高松矿田为例

刘春学, 倪春中, 燕永锋, 谭 喨

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

利用地质统计学方法模拟岩石裂隙网络的三维空间分布
——以云南个旧高松矿田为例

刘春学1), 倪春中2), 燕永锋2), 谭 喨1)

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

裂隙在地学的诸多领域中均具有重要意义, 其空间分布可以使用地质统计学方法进行模拟, 同时考虑裂隙的方向(走向和倾角)。利用序贯高斯模拟方法可以估计裂隙密度的空间分布, 并根据裂隙密度数值随机产生裂隙位置的空间分布。裂隙方向被划分成若干(非)均等的方向组, 将裂隙方向归属到其所属方向组,表示为由若干二值变量组成的指示形式, 0和1分别代表该裂隙方向不属于和属于该组。为了便于计算, 减少方向指示变量的成分数目, 使用主成分分析法求出方向指示变量的主成分, 用普通克里格法估计各主成分的空间分布。把估计结果反演为原始的指示形式, 并找出其中数值最大的方向组且将其赋值为1。按照对应方向组内裂隙方向的累积密度函数, 随机产生裂隙的方向。根据估计结果, 将符合一定距离和角度标准的裂隙元连接为一个裂隙面, 从而形成裂隙网络。根据在云南个旧锡矿高松矿田白云岩中进行裂隙网络模拟的应用, 可见该方法由于组合了序贯高斯模拟法、普通克里格法和主成分分析法, 可以较好地对岩石裂隙位置和方向进行合理的模拟。

裂隙网络; 三维空间分布; 地质统计学; 个旧锡矿

岩石中的裂隙网络分布在地学的诸多领域中都有重要的意义, 严重影响着建筑中岩石的稳定和油气、地下水、地热等地下资源的渗透和存贮(汪洋等, 2001; 曾普胜等, 2004; 张秋文等, 2004; 付广等, 2007; 田洪水等, 2007)。虽然近年来在不同的领域中都有许多方法用于模拟裂隙的分布(Snow, 1969; Rouleau et al., 1985; Long et al., 1987; Chilès, 1988; Acuna et al., 1995; Clemo et al., 1997; Riley, 2004),但是这些方法多数局限于用随机函数(高斯或泊松)生成裂隙面, 而对裂隙面的位置和方向考虑较少,与实际裂隙网络之间的差异较大。本文使用地质统计学方法模拟裂隙的三维分布, 同时考虑裂隙的方向(走向和倾向), 并应用于个旧锡矿高松矿田的裂隙网络模拟。

1 地质统计学模拟裂隙网络分布的方法

岩石裂隙本身具有多重属性, 诸如位置、方向、宽度、形状等; 不同的裂隙可能通过同一个空间位置; 裂隙网络经常表现出继承性和某种程度的尺度不变性。为了考虑裂隙的这些特点, 模拟裂隙的空间位置和方向、连接裂隙元, 可以使用数学地质方法(Koike et al., 2001, 2012)。该方法由普通克里格(OK: Ordinary Kriging)(Journel et al., 1978)、序贯高斯模拟(SGS: Sequential Gaussian Simulation)、主成分分析(PCA: Principal Component Analysis)等方法构成。

1.1 裂隙位置模拟

若将裂隙面定义为空间一定直径的圆盘, 则裂隙位置可以用其在三维空间的中心点来刻画, 裂隙密度(FD: Fracture Density)可以用单位体积、面积、长度内裂隙面、迹线、交切点的数目来表示。裂隙的空间分布可以根据裂隙密度数值用蒙特卡洛(Monte Carlo)法进行模拟。

虽然裂隙位置经常表现出空间偏倚和丛聚性,但是裂隙密度在多数情况下呈正态分布, 因此裂隙密度是模拟裂隙空间位置的有用工具。利用半变异函数阐明裂隙密度的空间关系, 进而在半变异函数模型的基础上应用OK或SGS估计裂隙密度的空间分布。在多数情况下, SGS在模拟裂隙密度的空间分布上表现出较好的效果(Koike et al., 1999, 2006; Liu et al., 2006), 这主要是因为SGS能减少克里格的平滑效果。

在进行估值之后, 研究区域内的每个待估单元格都被赋予一个裂隙密度值, 即该单元格内的裂隙数目。从而研究区域内裂隙位置的空间分布可以根据裂隙密度值应用蒙特卡洛法按照均匀分布随机产生。

1.2 裂隙方向模拟

由于构造应力的作用, 裂隙一般具有优势方向,因此裂隙方向不能简单的应用普通克里格法进行估计, 而需要应用主成分分析法进行估计, 可以按下列步骤进行。

将裂隙方向范围均等或不均等的分为n个互斥的方向组, 并根据裂隙方向所属的组将其表示为指示形式(g1g2… gn)。其中只有裂隙方向角度所在的组被赋值1, 代表该裂隙方向出现在该方向组内, 而其他组被赋值0。例如, 裂隙的走向范围[0, π]可以被分为4个方向组(EW NE NS NW), 其中EW代表[0, π/8)∪(7π/8, π), NE、NS、NW分别代表(π/8, 3π/8]、(3π/8, 5π/8]、(5π/8, 7π/8]。因此走向为π/4的裂隙被归入NE组, 对应的指示值为(0 1 0 0)。裂隙方向的组数可以自由决定, 但一般说来n=4就足够了; 如果考虑裂隙的空间方位, 即走向和倾向, 一般n=8是可以接受的。

由于裂隙方向的指示值包括多个成分, 主成分分析被用来减少成分的数目。一般情况下, 三个主成分即可在高置信度下代表所有成分变量的信息量,而三个成分变量对于计算机程序计算来说是可行的。

通过计算各个成分的半变异函数, 找出其中隐藏的空间变异规律, 进而用普通克里格法或序贯高斯模拟法估计各主成分的空间分布。从而在研究区域内的每个待估裂隙点上, 估计裂隙方向的所有主成分。

将每个位置上估计的裂隙方向的主成分值反演为指示形式。将具有最大值的成分用1来表示, 代表裂隙出现在该方向组的可能性最大; 其他成分赋值为0, 代表裂隙不在这些方向组中出现。例如, 在某一待估位置上的反演指示值为(0.5 0.3 0.1 0.9),将其中的最大值0.9替换为1, 其他的替换为0, 变成(0 0 0 1), 据此可以认为该裂隙方向最可能的方向组为NW方向组(5π/8, 7π/8)。

根据该方向组内实验累积分布函数(CDF: Cumulative Distribution Function), 在该方向组角度范围内随机产生裂隙的方向。

1.3 裂隙连接

对于估计的裂隙元, 根据给定的距离和方向角度标准, 可以将裂隙空间位置和方向相近的裂隙元连接起来, 视为属于同一个裂隙面, 按照一定的规则连接为一个裂隙面。真实的裂隙面本身是个复杂的空间面, 在不同的位置具有不同的方向, 因此可以将裂隙元连接为一个三维三角网; 也可以用裂隙元的平均走向和倾角连接为一个三维空间平面, 该平面的边界置于裂隙元在该平面上的投影位置。

如果某一裂隙元未能与其他裂隙元连接, 是孤立的, 可以根据其走向和倾角将之表示为一个直径为LL的空间圆盘。

2 个旧锡矿高松矿田应用

作为应用, 上述方法被用于模拟个旧锡矿高松矿田岩石裂隙网络的空间分布, 该研究区位于云南省东南部的个旧市附近(图1)。

2.1 高松矿田构造概况

高松矿田是云锡公司的主要生产矿山之一, 位于个旧锡矿东矿区北部, 处于矿区的五子山复背斜北段, 夹持于南北向个旧断裂、甲界山断裂与东西向个松断裂、背阴山断裂之间, 地质条件复杂。地表无矿体出露, 仅在芦塘坝、马吃水、高峰山三个地段的深部找到锡多金属矿床, 总体沿NE向芦塘坝断裂呈带状分布。高松矿区主要出露中三叠统地层, 以白云岩、白云质灰岩为主, 另有第三系泥岩、第四系残坡积物零星出露。

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

图1 个旧矿区高松矿田构造地质略图(据庄永秋, 1996, 修改)Fig. 1 Schematic geological map of the Gaosong orefield in the Gejiu tin mine (modified after ZHUANG, 1996)

构造是影响个旧锡矿高松矿田的矿体产出的重要因素之一。五子山复式背斜及其次级褶皱与北东向、近东西向断裂相互配置及其与花岗岩、有利地层的交割关系等, 不但为深部岩浆侵位提供了有利空间及成矿作用集中的场所, 并对矿田、矿床以至矿体起到了具体的定位作用(於崇文等, 1990; 庄永秋等, 1996; 刘春学等, 2003)。

2.2 样本裂隙数据

在高松矿田芦塘坝附近约EW 4 km× NS 5 km×H 0.8 km的范围内(图2), 根据总长度约5463 m的7条矿山生产坑道编录共搜集整理了12212条节理、裂隙、裂隙带等样本裂隙数据, 分布在1360、1540、1600、1720、1920等中段上。坑道主要呈北东和北西方向水平展布, 只有2条斜井,因此样本裂隙的产状多为垂直和陡倾斜、较少水平和缓倾斜(图2)。虽然裂隙的方向偏差可以在一定假设前提下进行校正, 但得到的结果差别不大, 因此直接利用样本裂隙数据模拟裂隙的三维空间分布。根据坑道和样本裂隙的空间分布区域, 限定研究区域的范围为4000 m(E-W)×5000 m (E-W)×600 m (Vertical), 待估单元格子大小设定为100 m× 100 m×25 m。

2.3 裂隙位置模拟

既然裂隙密度由单位长度内的裂隙数目决定,那么长度单位在计算裂隙密度时就非常重要, 并且影响计算得到的裂隙密度是否适用于普通克里格模型。实际上, 分别以5 m和10 m计算了裂隙密度, 从其频率分布图来看, 两种裂隙密度均服从对数正态分布。但从其半变异函数图和交叉验证系数来看,以10 m为单位的裂隙密度的半变异函数具有较大的变程和拱高, 得到的交叉验证系数也较大。因此,选取以10 m为单位计算裂隙密度(图3), 计算了其在NE、NS和垂直三个方向上的半变异函数, 同时计算了其在无方向情况下的半变异函数(图3)。从半变异函数的形状及其交叉验证系数来看, 无方向半变异函数可以提供较好的估计精度, 这也说明裂隙密度可以认为是各向同性的。以裂隙密度的各向同性半变异函数模型为基础, 使用SGS方法估计了裂隙密度的空间分布(图3)。在每个待估格子上, 利用随机函数产生与估计的裂隙密度数值一样的裂隙个数。

2.4 裂隙方向模拟

裂隙方向由走向(α)和倾角(β)构成, 表示为(α, β), 其中走向的范围为0到π, 倾角的范围为0到π/2。从个旧锡矿高松矿田搜集到的裂隙样本来看,裂隙走向主要为北西向, 北东向次之; 倾角较大,多在π/4以上(图4)。

首先裂隙方向在走向(0～2π)和倾角(0～π/2)范围内被分为16个相等的方向组。但由于各方向组内的裂隙数目太少, 很难计算半变异函数, 因此为了使各方向组内都能有合理的裂隙数目, 将裂隙方向分为了8组。

图2 样本裂隙位置俯视图(A)和立体图(B)Fig. 2 Bird view map (A) and stererogram (B) of sample fracture locations

图3 样本裂隙密度频率图(A)、半变异函数图(B)及估计裂隙密度空间分布图(C)Fig. 3 Histogram (A), semivariogram (B) of sample FD, and spatial distribution (C) of simulated FD

图4 样本裂隙方向频率分布图(A)及极点图(B)Fig. 4 Histogram (A) and diagram (B) of sample fracture directions

图5 样本裂隙方向主成分1(A)、2(B)、3(C)的半变异函数图Fig. 5 Semivariograms of the principal component 1(A), 2(B) and 3(C) from sample fracture directions

对应8个方向组, 每个裂隙方向均被转换为由8个二值数字(0和1)组成的指示形式, 例如样本裂隙方向(π/8, π/4)被转换成(1 0 0 0 0 0 0 0)。利用PCA分析指示变量, 求出指示值的主成分, 进而计算各个主成分的实验半变异函数, 并用球状半变异函数模型模拟出理论半变异函数(图5)。根据模拟得到的理论半变异函数, 利用普通克里格方法计算各个主成分的空间分布。将待估点上估计的主成分值反演为包含8个数值的指示形式, 其中最大值对应的方向组作为该点裂隙方向所属的组, 根据该方向组内样本裂隙走向和倾角的实验CDF, 用随机函数产生该裂隙的方向。

图6 裂隙元连接示意图Fig. 6 Schematic diagram showing connection of fracture elements

图7 个旧锡矿高松矿田裂隙网络空间分布俯视图(A)和立体图(B)Fig. 7 Bird view map (A) and stereogram (B) of simulated fracture networks in the Gaosong orefield of the Gejiu tin mine

2.5 裂隙元连接

根据上述方法估计得到的裂隙位置和裂隙方向,可以确定基本的裂隙元, 进而可以将空间分布接近的裂隙元连接为一个裂隙面。用空间面表示裂隙,可以根据下列距离和角度的标准将相近的裂隙元连接起来:

Dis≤PD; Dis为两个裂隙面之间的距离, PD为允许长度。

φ≤PA; φ为两个裂隙面之间的夹角, PA为允许角度。

将连接的裂隙投影到由其平均走向和倾角决定的平面中, 并将裂隙元的投影点连接为多边形(图6)。对于孤立的裂隙根据其走向和倾角表示成直径为10 m的圆盘。

图7显示了那些包含了60个裂隙元的裂隙面,不同的裂隙面用不同颜色表示。

2.6 结果及讨论

根据裂隙网络模拟及连接的结果, 可以看出,大裂隙面(包含60个以上裂隙元)的空间分布与地表观察到的大断裂之间具有较好的吻合关系, 体现了整个研究区的断裂分布特征, 可以较明显地看出断裂走向的变化。

模拟得到的大裂隙面主要集中研究区的中部和北东部, 多数呈北西向展布, 与研究区北西向断裂发育的情况一致; 东西向的裂隙面较少, 主要出现在研究区的南部, 与背阴山断裂比较吻合; 南北向的裂隙面在研究区的中部和北东部出现, 需要进一步分析研究区南北向的实测断裂; 另外模拟结果中还有倾角较缓的裂隙面出现, 值得开展进一步实地考察工作。但由于样本裂隙数据的采样偏差, 北东向的裂隙面很少, 只有在研究区的东南部出现了产状平缓的北东向裂隙面。

3 结论

利用SGS模拟裂隙位置的空间分布、利用主成分分析和普通克里格法模拟裂隙方向的空间分布,可以较适合地模拟裂隙的空间分布, 且能同时考虑裂隙的方向。从在个旧锡矿高松矿田白云岩中的实际应用结果来看, 模拟裂隙与实际裂隙之间无论在空间位置上, 还是在方向的频率分布上都具有较好的对应关系。

SGS方法比较适合模拟裂隙的位置, 可以较好地结合普通克里格的趋势性和随机过程的变化性。PCA可以合理地模拟裂隙方向的空间分布, 可以较好地体现裂隙方向的非均匀性, 从模拟裂隙显示的主要方向来看, 与实际断裂的空间展布方向一致。

付广, 吕延防, 杨一鸣. 2007. 我国大中型气田断裂输导天然气能力综合评价[J]. 地球学报, 28(4): 382-388.

刘春学, 秦德先, 党玉涛, 谈树成. 2003. 个旧锡矿高松矿田综合信息矿产预测[J]. 地球科学进展, 18(6): 921-927.

田洪水, 李洪奎, 王金光, 郭广军, 张增奇. 2007. 沂沭断裂带及其近区的地震成因岩石新认识[J]. 地球学报, 28(5): 496-505.

汪洋, 旸汪集, 熊亮萍, 邓晋福. 2001. 中国大陆主要地质构造单元岩石圈地热特征[J]. 地球学报, 22(1): 17-22.

於崇文, 蒋耀淞. 1990. 云南个旧成矿区锡石-硫化物矿床原生金属分带形成的动力学机制[J]. 地质学报, 64(3): 226-237.

曾普胜, 王海平, 莫宣学, 喻学惠, 李文昌, 李体刚, 李红, 杨朝志. 2004. 中甸岛弧带构造格架及斑岩铜矿前景[J]. 地球学报, 25(5): 535-540.

张秋文, 张培震, 王乘, MICHAEL A E. 2004. 断层间相互作用的触震与缓震效应定量评价[J]. 地球学报, 25(4): 483-488.

庄永秋, 王任重, 郑树培, 尹金明. 1996. 云南个旧锡铜多金属矿床[M]. 北京: 地震出版社.

References:

ACUNA J, YORTSOS Y. 1995. Application of fractal geometry to the study of networks of fractures and their pressure transient[J]. Water Resources Research, 31(3): 527-540.

CHILÈS J P. 1988. Fractal and geostatistical methods for modeling of a fracture network[J]. Mathematical Geology, 20(6): 631-654.

CLEMO T, SMITH L. 1997. A hierarchical model for solute transport in fractured media[J]. Water Resources Research, 33(8): 1763-1783.

FU Guang, LÜ Yan-fang, YANG Yi-ming. 2007. A Comprehensive Evaluation of Gas-transporting Capacity through Faults of Large and Medium Gas Fields in China[J]. Acta Geologica Sinica, 28(4): 382-388(in Chinese with English abstract).

JOURNEL A G, HUIJBREGTS CH J. 1978. Mining Geostatistics[M]. Salt Lake City: Academic Press.

KOIKE K, ICHIKAWA Y. 2006. Spatial correlation structures of fracture systems for identifying a scaling law and modeling fracture distributions[J]. Computers & Geosciences, 32(8): 1079-1095.

KOIKE K, KANEKO K. 1999. Characterization and modeling of fracture distribution in rock mass using fractal theory[J]. Geothermal Science & Technology, 6: 43-62.

KOIKE K, KOMORIDA K, ICHIKAWA Y. 2001. Fracture-distribution modeling in rock mass using borehole data and geostatistical simulation[M]. Proceedings of International Association for Mathematical Geology Cancun 2001 Conference, CD-Rom Press.

KOIKE K, LIU Chun-xue, SANGA T. 2012. Incorporation of fracture directions into 3D geostatistical methods for a rock fracture system[J]. Environmental Earth Sciences, 66(8): 1403-1414.

LIU Chun-xue, QIN De-xian, DANG Yu-tao, TAN Shu-cheng. 2003. Synthesis Information Based Mineral Resource Prediction of Gaosong Field in Gejiu Tin-deposit[J]. Advance in Earth Sciences, 18(6): 921-927(in Chinese with English abstract).

LIU C, KOIKE K, SANGA T. 2006. Geostatistical simulation of rock fractures distribution by considering directional elements[C]. Proceedings of International Association for Mathematical Geology 2006 Conference (Liège), CD-Rom Press.

LONG J C S, BILLAUX D M. 1987. From field data to fracture network modeling: an example incorporating spatial structure[J]. Water Resources Research, 23(7): 1201-1216.

RILEY M S. 2004. An algorithm for generating rock fracture patterns: mathematical analysis[J]. Mathematical Geology, 36(6): 683-702.

ROULEAU A, GALE J E. 1985. Statistical characterization of the fracture system in the Stripa granite, Sweden[J]. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts, 23(3): 353-367.

SNOW D T. 1969. Anisotropic permeability of fractured media[J]. Water Resources Research, 5: 1273-1289.

TIAN Hong-shui, LI Hong-kui, WANG Jin-guang, GUO Guang-jun, ZHANG Zeng-qi. 2007. New Recognition of Seismogenic Rocks in the Yishu Fault Zone and Its Periphery[J]. Acta Geoscientica Sinica, 28(5): 496-505(in Chinese with English abstract).

WANG Yang, WANG Ji-yang, XIONG Liang-ping, DENG Jin-fu. 2001. Lithospheric Geothermics of Major Geotectonic Units in China Mainland[J]. Acta Geoscientica Sinica, 22(1): 17-22(in Chinese with English abstract).

YU Chong-wen, JIANG Yao-song. 1990. The Dynamic Mechanisms of Primary Metal-zoning of Cassiterite-sulfids Deposits in The Gejiu Ore District, Yunnan Province[J]. Acta Geoscientica Sinica, 64(3): 226-237(in Chinese with English abstract).

ZENG Pu-sheng, WNAG Hai-ping, MO Xuan-xue, YU Xue-hui, LI Wen-chang, LI Ti-gang, LI Hong, YANG Zhao-zhi. 2004. Tectonic Setting and Prospects of Porphyry Copper Deposits in Zhongdian Island Arc Belt[J]. Acta Geoscientica Sinica, 25(5): 535-540(in Chinese with English abstract).

ZHANG Qiu-wen, ZHANG Pei-zhen, WANG Cheng, MICHEAL A E. 2004. Interaction of Active Faults and Its Effect on Earthquake Triggering and Delaying[J]. Acta Geoscientica Sinica, 25(4): 483-488(in Chinese with English abstract).

ZHUANG Yong-qiu, WANG Ren-zhong, ZHENG Shu-pei, YIN Jin-ming. 1996. Tin-copper-polymetallic Deposit in Gejiu, Yunan[M]. Beijing: Seismological Press(in Chinese).

Three-dimensional Simulation of Rock Fractures by Geostatistical Method: A Case Study of Gaosong Field in Yunnan Province

LIU Chun-xue1), NI Chun-zhong2), YAN Yong-feng2), TAN Liang1)
1) School of Urban and Environment, Yunnan University of Finance and Economics, Kunming, Yunnan 650221; 2) Faculty of Land Resources Engineering, Kunming University of Science and Technology, Kunming, Yunnan 650093

The simulation of fracture distribution is an important problem in various fields of geosciences. To simulate the spatial distribution of fracture networks, this paper proposes a geostatistical method in consideration of their directions (strikes and dips). Fracture locations are generated randomly, following fracture density values assigned by sequential Gaussian simulation method. Fracture direction is divided into n equal (or unequal) groups, and sample fracture directions are assigned to its corresponding group. Then sample fracture directions are transformed into indicators consisting of n binary variables, where 1 and 0 represents belonging to and not belonging to this group. For calculation convenience, the indicator number is reduced by using the principal component analysis. Then ordinary kriging is employed to estimate the distributions of these principal components. The results are inversed to the original indicator form, and the biggest one is assigned as 1 while the others areassigned as 0. Fracture directions are generated randomly by using the cumulative distribution function of the biggest group. Based on these simulated results, fracture elements can be determined with location and direction. At last, fracture elements within the angle and distance tolerances are connected to be one fracture. The case study of the fracture data in Gejiu dolomite of southeast Yunnan Province shows that the combination of sequential Gaussian simulation, ordinary kriging and principal component analysis can provide a reasonable simulation result for locations and directions of fractures.

fracture network; three dimensional distribution; geostatistics; Gejiu tin mine

O242.1; P628.2

A

10.3975/cagsb.2013.03.10

本文由国家自然科学基金项目“地学中方向性变量的多尺度空间分布模拟”(编号: 40902058)资助。

2012-12-05; 改回日期: 2013-02-25。责任编辑: 魏乐军。

刘春学, 男, 1975年生。博士, 教授。主要从事裂隙模拟和水环境研究。通讯地址: 650221, 云南省昆明市龙泉路237号云南财经大学城市与环境学院。电话: 0871-5183357。E-mail: chunxueliu@hotmail.com。