裂隙网络的三维随机模拟方法研究

沈瑞文

(安徽理工大学 地球与环境学院，安徽 淮南 232001)

裂隙网络的模拟研究在地质学领域的诸多方面都有十分重要的应用，裂隙网络的准确刻画和模拟是地质学领域具有挑战性的问题之一[1]。利用随机方法模拟裂隙网络最早于20世纪80年代末由Neuman提出，该方法早期多应用于模等效连续介质模型，90年代后期逐步应用于裂隙网络的随机模拟[2]。目前，裂隙网络模拟技术已逐渐成熟，其成果已在工程上得到了广泛的应用。陶凯等利用加权Voronoi图法生成二维裂隙网络模型，并应用在膨胀土的力学性质和渗流特性[3]；王晋丽等利用Monte-Carlo随机模拟技术生成二维裂隙网络渗流模型，研究了裂隙网络中水流在不同边界条件下的流动特征[4]；敖雪菲等在研究坝基裂隙岩体灌浆时，基于Monte-Carlo方法建立了裂隙岩体三维网络模型，结果表明该法能有效反应裂隙真实特征并获得贴近实际的灌浆结果[5]；张弛等利用布尔模拟和多点统计学方法进行了裂隙网络的三维模拟，并对比两种模型的精度[6]；李瑞金等基于Monte-Carlo方法实现裂隙的三维随机模拟，并应用于坝基多孔分序灌浆数值模拟[7]。

近年来，随着裂隙网络模拟研究在采矿设计以及地下工程施工中应用的增加，对裂隙网络模拟的精度要求越来越高。以往根据裂隙的定量数据对裂隙进行二维模拟的方法仅仅对模拟对象进行平面展布，没有考虑模拟对象方位角度的变化，与实际裂隙网络之间的差异较大，难以达到高精度地质模型的要求[8]。三维地质随机模拟技术相较于二维模拟具有诸多优点，它在原二维模拟的基础上考虑了裂隙的位置和方向(走向和倾角)，有利于对地质变量各属性间的关系进行充分研究，确定裂隙的变化趋势，建立更加符合实际情况的裂隙网络，为采矿设计、地质资源预测、地下工程施工、地质灾害预报等领域的研究和工作提供科学的定量依据[9]。本文在裂隙变量数据的基础上，采用Monte-Carlo和地质统计学相结合的三维随机模拟方法模拟出裂隙元，并同时考虑裂隙的距离和角度，按照一定的标准连接成裂隙面，形成高精度的、与实际裂隙吻合程度较高的裂隙网络。

1 三维随机模拟原理

近年来，在不同领域有许多方法用来模拟裂隙网络的分布，常见的方法有加权Voronoi图法、序贯指示模拟、布尔模拟、多点统计学模拟、拉丁超立方抽样模拟(LHS)、Monte-Carlo等等[10]。Monte-Carlo模拟是通过一定的随机数生成方法，生成服从某一随机变量的概率分布形式的随机数序列[11]。裂隙具有方向、位置、形状等多重属性，它在岩体中的分布是随机的、复杂的[12]，要对裂隙网络进行模拟具有一定的难度，但是大量实践证明裂隙的定量数据(密度)往往服从一定的概率分布，因此，可以运用Monte-Carlo来研究裂隙分布特征以及模拟裂隙网络分布[13]。地质统计学的基本分析方法是克里金法，克里金是一种确定性内插方法，通过它计算得到的实现期望偏差相对最小。克里金方法是建立在随机地质介质模型基础之上(即变差函数基础上)的确定性求解，可以求出与未知的客观实际值之间偏差最小的唯一解。将Monte-Carlo和地质统计学结合在一起，流程如图1所示，对裂隙网络进行三维模拟研究，得出的结果与实际裂隙之间的吻合性较高，同时可以更好地为后续的采矿设计以及地下工程施工提供精度较高的依据[14]。

图1 裂隙网络模拟流程

2 裂隙网络非定量数据的模拟方法

根据上文所述的原理，归纳出裂隙网络模拟实现的步骤：

1)裂隙位置模拟。在裂隙网络的三维随机模拟过程中，裂隙密度的计算占有首要位置，它反映了裂隙在研究区域中发育的密集程度,密度是定量数据，以往通常采用线密度或面密度等二维形式进行计算,这里对裂隙密度用三维形式——体密度进行计算，体密度是指研究区域内单位体积分布的裂隙条数[15]，其计算公式为

(1)

式中：n为裂隙总数目；ri为裂隙平均直径；V为取样体积。

对于研究区域内的每一个空间单元格，可以计算落入其中的裂隙中心点的数目，从而根据单元的大小得到各个单元格的裂隙密度。

Monte-Carlo随机模拟方法的基本原理是利用[0，1]区间标准均分随机数，根据由已知裂隙几何参数的密度分布函数求得的抽样公式来获得服从给定裂隙几何参数分布形式的随机变量的模拟模型，其抽样方法中直接抽样方法最为常用和有效，具体实施方案：设随机变量x服从分布函数y(x)和累积分布函数Y(x)，Y(x)的值域为[0，1]，ti作为Y(xi)的函数值，则x和t之间的对应关系为[16]

(2)

其反函数为

x=Y-1(t).

(3)

通过计算机程序产生一系列均匀分布随机数x1，x2，…，xn，将这些均匀分布随机数代入式(2)，即得到随机变量X的模拟模型[17]

Xi=Y-1(ti)，i=1,2,…,n.

(4)

式中：n为试验次数，实践表明，试验次数越多，X的频率分布越接近其真实的概率分布，在实际中n一般取频率分布收敛时所对应的试验次数，在三维随机模拟研究过程中，当模拟的次数足够多时，就可以得到相对比较精确的概率值[18]。

裂隙位置的模拟是根据已经计算出的具有正态分布特征的密度值，运用Monte-Carlo和地质统计学相结合的方法，在密度的半变异函数模型的基础上，应用普通克里格对裂隙密度的空间展布进行模拟，并随机模拟出单位体积内的裂隙条数，最后得到裂隙位置空间分布图。

2)裂隙方向模拟。走向和倾角在裂隙几何特征中占有重要位置。裂隙在构造应力的作用下会产生优势方向，根据走向玫瑰花图、倾角分布的直方图对其产状分布特征进行分析，用n个互斥的方向组来表示裂隙的方向，裂隙方向用(t1t2…tn)形式表达。若裂隙方向在该方向组内出现，则这一组的值为1，同时，其他组因裂隙方向不在这些组内，则赋值为0[19]，从而地质表示的裂隙方向即可以实现向指示形式表示的转换。

主成分分析实质上是一种数学变换的方法，它把一组相关变量通过线性变换转换成另一组不相关的变量，将这些新的不相关的变量按照方差依次递减的顺序排列，就形成所谓的主成分，使得第一主成分具有最大的方差，第二主成分的方差次之，并且和第一主成分不相关，依次类推，形成p个主成分[20]。

运用主成分分析法对裂隙的指示值进行计算，得出裂隙方向组的主成分，并计算主成分的实验变异函数，再用普通克里格法模拟出计算得到的主成分的空间分布，用n个0和1组合成的指示形式将预估点处的主成分值反演出来，其中，最大值对应的方向组作为该点裂隙方向所属的组, 根据该方向组内样本裂隙走向和倾角的累积分布函数, 用随机函数产生该裂隙的方向[21]，如图2所示。

图2 裂隙方向模拟流程

3)裂隙面的建立。对已经模拟出的裂隙元，以一定的距离和角度方向为标准，视空间分布位置、方向相同或相近的裂隙同属于一个相同的裂隙面，将其以一定的准则连接为裂隙面[22]。

令裂隙中心点之间的距离为Lc，裂隙与其中心点连线的夹角为α1，α2，如图3所示，设L为裂隙元连接的最大距离，一般规定为2～3个单元格大小；α为裂隙元连接的最大角度，一般规定为5°～10°。将符合式(5)的裂隙连接成同一个裂隙面[23]

图3 裂隙元连接

Lc

(5)

3 研究实例

由于朱仙庄矿区煤层主要位于二叠系地层内部，本文进行裂隙模拟的地层即为朱仙庄煤系地层——二叠系。对朱仙庄矿区内的实际裂隙数据进行搜集，对搜集到的样本数据进行分析，以10 m为单位，统计单元格内的裂隙数目，并以此计算出裂隙密度，如图4所示,得出的数据符合正态分布，并根据实际裂隙的展布情况，做出走向玫瑰花图，如图5所示，可以看出裂隙主要分布在NW、NE等方向。

图4 裂隙密度频率

图5 走向玫瑰花图

根据裂隙密度，将Monte-Carlo模拟和地质统计学方法相结合，在密度的半变异函数模型基础上，应用普通克里格模拟出裂隙空间分布，用普通克里格对裂隙密度的空间展布进行模拟，并随机模拟出单位体积内的裂隙条数，最后得到裂隙位置空间分布图，如图6所示。

图6 裂隙密度的半变异函数(A)、裂隙密度空间分布(B)和裂隙位置空间分布(C)

将裂隙方向在走向a∈[0°,360°]和倾角b∈[0°,360°]的范围内分为6个方向组：

T1{a∈[0°,60°),b∈(0°,90°]},

T2{a∈[60°,120°),b∈(0°,90°]},

T3{a∈[120°,180°),b∈(0°,90°]},

T4{a∈[180°,240°),b∈(0°,90°],

T5{a∈[240°,300°),b∈(0°,90°]},

T6{a∈[300°,360°),b∈(0°,90°]}.

对这6个方向组进行处理，将其用二值数字0和1组成的指示形式表示出来。利用主成分分析这6个指示变量，求出3个主成分(PC1，PC2，PC3)，对其进行计算并分别模拟出半变异函数模型，如图7所示。基于半变异函数模型应用普通克里格模拟出PC1，PC2，PC3的空间分布，如图8所示。按照6个0和1构成的变量反算出待估裂隙的主成分值，则该裂隙方向即为最大值对应的方向组，根据模拟得到的裂隙方向组，应用Monte-Carlo和地质统计学相结合，模拟出裂隙方向。

图8 3个主成分对应的空间分布(Ⅰ、 Ⅱ、 Ⅲ)

图7 3个主成分对应的半变异函数(Ⅰ、 Ⅱ、 Ⅲ)

根据三维随机模拟得到的裂隙元，在综合考虑其距离和角度的基础上，将属于同一裂隙面的裂隙元连接起来，形成裂隙面，如图9所示。可以看出，研究区域内的大裂隙面主要分布在中部、东南部，展布较为均匀，裂隙的发育方向主要在NW和NE。在应用Monte-Carlo进行模拟时，模拟次数越多结果越接近实际情况，而在实际模拟过程中，模拟次数不可能趋于无穷多，同时在收集实际裂隙数据时存在一定的误差，因此，模拟结果与实际情况具有一定的偏差。综合考虑多次模拟结果，得到与实际裂隙的走向趋势和空间分布具有较好的吻合性的模型。在模拟出的裂隙三维网络的基础上，可以清晰地看出研究区域内断裂走向变化的趋势，从而为进一步的采矿设计、地下工程施工等提供依据。

图9 实际裂隙分布(A)和模拟的裂隙网络空间分布(B)

4 结束语

运用Monte-Carlo和地质统计学相结合的方法对裂隙网络进行三维随机模拟，在Monte-Carlo模拟方法的实现过程中，会产生多个等概率的反映裂隙网络空间结构的模拟结果，各个结果间的差别反映了模拟次数、累积概率以及资料不足等原因而引起的不确定性， 综合考虑这些模拟结果，可以得到与实际裂隙展布情况较一致的裂隙网络。在对朱仙庄含煤系地层的裂隙网络模拟研究上，可以清楚地看出模拟出的裂隙与原实际裂隙在空间位置分布和方向变化趋势上都具有较好的吻合性，进一步说明在运用Monte-Carlo模拟方法的基础上，结合地质统计学，综合考虑裂隙的方向，有利于实现高精度裂隙网络模型的模拟，对合理评价矿区内裂隙结构特征具有重要意义，同时也能够为采矿设计和地下工程施工提供精度较高的依据。