随机分布三维裂隙网络中裂隙对相交判断方法

摘要

随机分布三维裂隙网络中裂隙对相交判断方法，包括：将三维空间中的每个裂隙用圆盘模拟，根据给定的随机分布裂隙几何参数生成随机分布的裂隙网络；用轴向包围盒代替裂隙圆盘；将整个裂隙网络所在的岩块分成小块，块内包含的裂隙圆盘轴向包围盒和在切面上相交的裂隙圆盘轴向包围盒归入该小块；每一小块作为独立单元，块内根据裂隙圆盘轴向包围盒的空间特征，采用空间数据索引算法构建空间裂隙索引以快速进行初步相交判断，得到块内所有可能相交的裂隙对。进一步包括对所有可能相交的裂隙对进行精确判断的步骤。此方法速度快、成本低，能够大大减少裂隙岩体渗流计算量，明显缩短裂隙岩体流体渗流和溶质运移数值模拟计算时间。

说明书

技术领域

本发明涉及随机分布三维裂隙网络中裂隙对的相交判断方法，具体涉及裂 隙对的初步与进一步精确相交判断方法。

背景技术

裂隙化岩体广泛存在于地表浅层，是主要的地下流体渗透介质之一。由不 同产状的裂隙构成的裂隙网络形成了地下水及溶质和地下油气资源在裂隙化岩 体中运移的主要甚至唯一通道。自然界中的裂隙错综复杂，长短不一，裂隙间 并非完全连通，总是存在着那些与其他裂隙面均不相交的孤立裂隙，因此判断 裂隙是否相交是各种渗流计算的基础。首先，自然界中裂隙面的数量巨大，当 裂隙面的个数到达几十万级别时，相应的离散裂隙网络渗流与污染物运移数值 计算量也相当大。通过渗流路径的准确识别，可获得对渗流与污染物运移有贡 献的裂隙，而去除所有无贡献的裂隙如孤立的裂隙，进而有效减少相应的数值 计算负荷，提高渗流等计算的效率和速度。通过渗流路径的初步识别，可获得 所有对渗流与污染物运移可能有贡献的裂隙，去除一定数量的绝对无贡献的裂 隙，进而也可减少相应的数值计算负荷，提高渗流计算的效率和速度。需要说 明的是，裂隙网络中裂隙对的初步与精确相交判断分别建立了相应的渗流路径 初步识别与精确识别的基础。第二，裂隙张量是裂隙的一个重要参数，有学者 认为可根据连接的裂隙张量分量确定方向岩块水力传导系数并估计典型单元体 (REV)值，而在计算裂隙张量时，也需要查明裂隙的构造，进行裂隙网络的 相交判断。因此，开发出高效、稳定的三维裂隙网络中裂隙对的相交判断算法 成为裂隙渗流计算至关重要的问题。

1985年Stauffer提出了逾渗理论(percolation theory)，指出了4类对水 流没有贡献的裂隙。研究人员(Wang EZ，Yang CT.，1992；Yu QC，Huang XX.， 1992；Lu B，Chen JP，Wang LK.，2002；Zhao HL，Chen JP.2005；Yu QC， Wu X，Ohnishi Yuzo.，2006)对二维或三维裂隙网络的流动通道识别方法进 行了研究，但均为方法的探讨，且没有在数据量较大的情况下对方法进行测试， 因此没有从应用角度发现和解决算法计算量大的问题。研究人员(LIU HuaMei and WANG MingYu，2010)根据图论的相关知识，将三维裂隙网络概化为非连 通图，采用计算机中图的数据结构表示三维裂隙网络数据，编制相应算法的计 算机程序进行了裂隙岩体沿水力梯度方向渗流路径的识别。发现在此算法中， 裂隙网络进行裂隙对相交判断的前期工作时间开销相当大，尤其是对于尺度较 大、裂隙数目较多的情况，需要较长运算时间和占用较多计算机资源。

从以上的分析中可以看出现有的三维随机分布裂隙网络中裂隙对相交判断 方法的明显不足，使用传统的方式进行大尺度、数量众多的三维随机分布裂隙 网络数据是不现实的，这极大的限制了裂隙岩体渗流研究的尺度和研究可行性， 很难推广到实际应用中。所以在设备资源和计算能力有限的条件下，如何实现 一个快速的、有效的、并且消耗资源少的三维随机分布裂隙网络中裂隙对相交 判断方法成为一个重要的问题。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种速度快、成本 低的随机分布三维裂隙网络中裂隙对相交判断方法。此方法能够大大提高随机 分布三维裂隙网络中裂隙对相交判断算法效率，从而大大减少裂隙岩体渗流计 算量，明显缩短裂隙岩体流体渗流和溶质运移数值模拟计算时间，使大尺度裂 隙网络渗流与溶质运移模拟研究与工程应用的可行性显著提高。

本发明的技术解决方案是：

随机分布三维裂隙网络中裂隙对相交判断方法，包括以下步骤：

将三维空间中的每个裂隙用圆盘模拟，根据给定的三维随机分布裂隙几何 参数生成随机分布的裂隙网络；

将裂隙圆盘的外接球体的外接正立方体作为该裂隙圆盘的轴向包围盒，用 轴向包围盒代替三维空间中的裂隙圆盘；

将整个裂隙网络所在的岩块在三维空间中分成小块，每一小块内包含的裂 隙圆盘轴向包围盒和在切面上相交的裂隙圆盘轴向包围盒归入该小块；

每一小块作为独立单元，块内根据裂隙圆盘轴向包围盒的空间特征，采用 空间数据索引算法构建空间裂隙索引以快速进行初步相交判断，得到块内所有 可能相交的裂隙对。

在得到块内所有可能相交的裂隙对之后，可以进一步包括以下步骤：

将每块所有可能相交的裂隙对进行对接合并，去除其中重复的裂隙对，对 剩余可能相交的裂隙对进行精确相交判断，得到整个随机分布三维裂隙网络中 全部相交的裂隙对；其中精确判断具体为：根据已知裂隙圆盘的几何参数，确 定裂隙圆盘所在的平面，判断空间两个圆盘是否相交。

在得到块内所有可能相交的裂隙对之后，还可以进一步包括以下步骤：

对每块内所有可能相交的裂隙对进行精确相交判断，得到每块内所有相交 的裂隙对，将各块所有相交的裂隙对进行对接合并，去除其中重复的裂隙对， 得到整个随机分布三维裂隙网络中全部相交的裂隙对；其中精确判断具体为： 根据已知裂隙圆盘的几何参数，确定裂隙圆盘所在的平面，判断空间两个圆盘 是否相交。

所述将三维空间中的每个裂隙用空间圆盘模拟，根据给定的三维随机分布 裂隙几何参数生成随机分布的裂隙网络的步骤；具体包括：将每个裂隙模拟为 一个空间薄圆盘，其几何参数包括三维空间中裂隙圆盘的圆心坐标、倾向、倾 角、裂隙圆盘半径、裂隙密度；采用蒙特卡罗方法在三维空间中产生多个具有 不同大小和产状的薄圆盘。(产状是地质术语，指一个地质结构平面的倾向与倾 角。)

所述将整个裂隙网络所在的岩块在三维空间中分成小块，每一小块内包含 的裂隙圆盘轴向包围盒和在切面上相交的裂隙圆盘轴向包围盒归入该小块的步 骤，具体为：用一个切面将研究区域分割成若干小块，以块内包含的裂隙圆盘 圆心数作为描述块的基本大小；

分块方法采用三方向分块法，按照切割方向在三维空间XYZ三方向分块， 每隔规定的尺寸分为一块；或者

分块方法采用单方向分块法，按照切割方向在三维空间XYZ单方向分块， 即在空间上沿着XYZ三个方向中的任意一个方向分割。

所述每一小块作为独立单元，块内根据裂隙圆盘轴向包围盒的空间特征， 采用空间数据索引算法构建空间裂隙索引以快速进行相交判断，得到可能相交 的裂隙对的步骤，包括空间对象的插入操作、空间对象的删除操作和空间对象 的搜索操作。

所述空间对象的插入操作，具体包括：

(1)为新的空间对象寻找一个合适的叶子结点：选择合适的叶子结点L存放 新的空间对象；

(2)将新的空间对象记录在叶子结点中：若L中尚有空位(即L中的单元个 数小于L结点记录数目的最大值M)，则在L中记录新的空间对象的最小边界矩 形MBR和结点编号ID；若L中没有空位，则将L结点中的单元连同新的空间 对象分裂成两部分，分别记录在结点L和新生成的结点LL中；

(3)调整树的结构：从叶子结点L开始，对空间裂隙圆盘索引进行调整；

(4)生成新的根结点：若子结点的分裂导致根结点中单元个数超过M，则将 根结点分裂成两个结点，并再生成一个新的结点作为空间裂隙圆盘索引的根结 点，原根结点分裂成的两个结点作为其子结点；此时空间裂隙圆盘形成的索引 增加一层。

所述步骤(1)为新的空间对象寻找一个合适的叶子结点：选择合适的叶子结 点L存放新的空间对象；进一步包括：

(1.1)初始化：令N为空间裂隙圆盘索引的根结点；

(1.2)判定该结点是否为叶子结点：若N是叶子结点，则返回N，N即为找 到的叶子结点；

(1.3)选择合适的子树：若N不是叶子结点，则将新的空间对象依次加入到 该结点的各单元，分别计算出该结点所在区域I的面积，其中面积增长最小的 单元所在的子树即为合适的子树；令N为该单元所指向的子结点，转到步骤(2) 继续。

所述步骤(3)调整树的结构：从叶子结点L开始，对空间裂隙圆盘索引进行 调整；进一步包括：

(3.1)初始化：令N为叶子结点L，NN为因新的空间对象插入而分裂生成 的新的结点；

(3.2)判断N是否为根结点：若N是根结点，则停止返回；

(3.3)调整父结点相应单元的I：调整结点N的父结点P中相应单元的I，使 其正好包含N结点中的所有单元；

(3.4)根据需要进一步分裂父结点：若存在分裂生成的新的结点NN，则在N 的父结点P中加入指向该结点的单元，若P结点中的单元的个数超过了M，则 需要分裂结点P，生成一个新的结点PP，并令N＝P，NN＝PP；转到步骤(3.2) 继续。

所述空间对象的删除操作，具体包括：

(1)确定删除空间对象所在的叶子结点：通过空间裂隙圆盘索引的搜索操 作，确定删除空间对象所在的叶子结点L；

(2)删除空间对象：将结点L中的空间对象所在的单元删除；

(3)调整空间裂隙圆盘索引：从叶子结点L开始调整空间裂隙圆盘索引的操 作；

(4)调整根结点：若空间裂隙圆盘索引调整后，根结点只有一个孩子，则将 根结点删除，让其孩子结点成为空间裂隙圆盘索引的根结点，此时空间裂隙圆 盘索引减小一层。

所述步骤(3)调整空间裂隙圆盘索引，具体包括：

①初始化：令N为被删除单元的叶子结点L，Q为结点的集合，用以存放被 删除的结点，初始化为空；

②寻找N在父结点中的相应的单元：如果N为根结点就转向步骤⑥；反之 查找其父节点P，找到与结点N相对应的单元EN；

③删除单元数小于m的结点：若结点N的单元数小于m，则在结点P中 删除单元E_N，并将结点N加入集合Q中；

④调整MBR：若结点N没有被删除，则调整单元E_N中的I，使其正好包 含结点N中的所有的单元；

⑤调整上一层：令N＝P，转到步骤②继续；

⑥被删除的结点重新插入到空间裂隙圆盘索引中：采用空间裂隙圆盘索引 的插入操作，将Q中的所有的被删除的结点重新的插入到空间裂隙圆盘索引中。

所述空间对象的搜索操作，具体为：

设搜索的区域为S，则搜索区域S内的空间裂隙的过程如下：

(1)子树的搜索：从空间裂隙圆盘索引的根结点T开始，如果T不是叶子结 点，那么依次判断该结点中各单元的I与区域S的空间位置关系，若I与搜索 区域S相交，则该单元所指向的结点为待搜索的根结点，转到本步骤开始处继 续子树的搜索；若T是叶子结点则转至步骤(2)；

(2)叶子结点的搜索：若T是叶子结点，则依次判断其中的空间对象与搜索 区域S之间的空间位置关系，如果空间位置对象落在搜索区域S以内，那么表 明其满足搜索条件。

本发明与现有技术相比具有如下优点：本发明综合多种交叉学科的方法， 在随机分布三维裂隙网络数据尺度较大、数量较多的情况下，通过简化问题以 及降低问题规模等多方面入手，采用空间圆盘外接轴向包围盒、沿着三维坐标 轴方向进行分块和构建空间索引等算法提出一种快速而有效的随机分布三维裂 隙网络中裂隙对初步相交判断方法，根据需要，可以此为基础，进一步完成裂 隙对的精确相交判断，这些为后续识别相应渗流路径做准备。初步相交判断得 到的结果不仅可以作为精确相交判断的基础，也可以直接用于不要求绝对精确 相交判断结果的渗流计算，二者皆大大减少了不必要的渗流计算负荷。该方法 运算速度快，不需要依赖昂贵的专业设备和复杂的后续处理，节约硬件成本。

附图说明

以下将结合附图对本发明的具体实施方式进行说明。

图1为根据本发明的随机分布三维裂隙网络中裂隙对相交判断方法流程 图。

图2为三维裂隙网络模拟示意图。

图3为三维裂隙网络相交判断算法的两种分块方法示意图，其中(a)为 三方向分块，(b)为单方向分块。

图4为空间裂隙圆盘的轴向包围盒分布及其形成的索引示意图。

图5为第1套和第7套裂隙系统三种方案的最佳时间对比，其中(a)为 第1套裂隙系统，(b)为第7套裂隙系统。

图6为不同分块方法及块内不同算法运行结果对比，其中(a)为第1套裂 隙系统，(b)为第2套裂隙系统，(c)为第3套裂隙系统，(d)为第4套裂隙 系统，(e)为第5套裂隙系统，(f)为第6套裂隙系统，(g)为第7套裂隙系 统，(h)为第8套裂隙系统，(i)为第9套裂隙系统。

具体实施方式

如图1所示，本发明提供的随机分布三维裂隙网络中裂隙对相交判断方法， 包括以下步骤：

步骤一：将三维空间中的每个裂隙用空间圆盘模拟，根据给定的三维随机 分布裂隙几何参数采用蒙特卡罗方法生成随机分布的裂隙网络。

在裂隙岩体中的裂隙网络建立的圆盘模型中，假定岩体中裂隙发育呈圆盘 状，每个裂隙被模拟为一个空间薄圆盘，其几何参数包括三维空间中裂隙圆盘 的圆心坐标x，y，z，倾角α，倾向β，裂隙圆盘半径R，裂隙密度ρ。三维裂隙 网络的模拟就是在三维空间中产生多个具有不同大小和方向的薄圆盘。

对于地质上为数众多的尺寸为几十厘米到十几米的裂隙，因为不可能逐个 测量，只能做到相对数量的采样测量，要得到其几何参数，可采用统计分析的 方法，具体步骤如下：

①根据统计分析理论得到裂隙的统计规律，每种属性的具体统计分布形式 和统计参数的分布类型可选，用概率分布来表征，所述概率分布类型为均匀分 布、负指数分布、正态分布、对数正态分布、伽马分布和Fisher分布中的一种；

②在已知裂隙网络各几何参数的概率分布类型后，对这些概率分布进行随 机抽样，得到单个裂隙的单个几何参数，将这些几何参数组合起来，构成具体 的裂隙网络。

根据已知概率分布的随机抽样，用蒙特卡罗方法模拟具体的裂隙，在研究 域内得到一个与实际裂隙化岩体具有统计意义上等价的具体的裂隙网络。

步骤二：用轴向包围盒代替三维空间中的裂隙圆盘，以达到简化问题的目 的。

计算三维空间中裂隙圆盘的外接球体的外接正立方体的最小坐标和最大坐 标，即最靠近坐标轴原点值和最远离坐标轴原点值，得到轴向包围盒的对角坐 标。以此作为后续判断的依据。

设A，B为两个裂隙圆盘的三维包围盒，X为投影轴，CA，CB分别为A， B的中心点，PA，PB分别为点CA，CB在X轴上的投影。RA，RB分别为 裂隙圆盘包围盒A，B在X上的投影。若RA+RB＜PAPB，则在轴向上A和B不 相交，反之在轴向上A和B邻接或相交。当裂隙圆盘包围盒A，B在三条坐标 轴向上的投影均相交时，则A，B相交。因此裂隙圆盘间的相交测试最多只需6 次比较运算。

步骤三：将整个裂隙网络所在的岩块在三维空间中分成小块，每一小块内 包含的裂隙圆盘的轴向包围盒和在切面上相交的裂隙圆盘的轴向包围盒归入该 小块。

分块的思路是将图的规模缩小，化整为零。分块的目的是减少单次处理的 裂隙数目，分批次处理以加快速度，即把无向非连通图在空间上人为分割，降 低n值，这样增加的时间复杂度为空间分割图的时间和识别路径后对图的拼接， 却大大降低原问题的复杂度。分块实际上是用一个切面将研究区域分割成若干 块，以块内包含的圆心数作为描述块的基本大小，而根据不同裂隙系统的构造， 实际计算的块内裂隙还要包含分块时与切面相交的那些裂隙。

具体的分块方法可按照切割方向在三维空间XYZ三方向分块，也可以单方 向分块，如图3所示。不同的分块方法对算法运行的结果会产生一定的影响。 本发明采用的分块方法有两种，一种是在空间上沿一个方向分割，每隔规定的 尺寸就分为一块，这样可以降低每块的规模。另一种分块方法是在空间上对XYZ 三个方向均可以分割。经过测试，三方向分块方法比单方向分块在切面上产生 的裂隙数目少，也就是说，在进行整个裂隙岩块的相交判断时，重复计算的裂 隙数目少，由此使得三方向分块在运算效率上的提高，但是当研究区域地形为 长条形时，单方向分块具有本质上的优势，因此本发明中所述一般通用的有效 方法为三方向分块，在特殊情况下，单方向分块方法也是一种有效方法。

步骤四：每一小块作为独立单元，块内根据裂隙圆盘轴向包围盒的空间特 征，采用空间数据索引算法构建空间裂隙索引以快速进行初步相交判断，得到 块内所有可能相交的裂隙对。

计算三维空间中裂隙圆盘的外接球体的外接正立方体的最小坐标和最大坐 标，即最靠近坐标轴原点值和最远离坐标轴原点值，得到轴向包围盒的对角坐 标。

设A，B为两个裂隙圆盘的三维包围盒，根据裂隙圆盘的几何参数，可以 得到轴向包围盒的8个顶点坐标以及对角坐标，X为投影轴，CA，CB分别为 A，B的中心点，PA，PB分别为点CA，CB在X轴上的投影。RA，RB分 别为裂隙圆盘包围盒A，B在X上的投影。若RA+RB＜PAPB，则在轴向上A 和B不相交，反之在轴向上A和B邻接或相交。当裂隙圆盘包围盒A，B在 三条坐标轴向上的投影均相交时，则A，B相交。因此裂隙圆盘间的相交测试 最多只需6次比较运算。

如图4所示，采用空间数据索引算法，构建三维裂隙圆盘索引以快速进行 相交判断得到可能相交的裂隙对的主要步骤包括：

1.插入算法：

空间裂隙圆盘索引的插入算法要能够保证树的平衡，使所有树的叶子结点 都处在树的同一层上。

(1)为新的空间对象，寻找一个合适的叶子结点：选择合适的叶子结点L来存 放新的空间对象。

(1.1)初始化：令N为空间裂隙圆盘索引的根结点；

(1.2)判定该结点是否为叶子结点：若N是叶子结点，则返回N，N即为找 到的叶子结点；

(1.3)选择合适的子树：若N不是叶子结点，则将新的空间对象依次加入到 该结点的各单元，分别计算出该结点所在区域I的面积，其中面积增长最小的 单元所在的子树即为合适的子树；令N为该单元所指向的子结点，转到步骤(2) 继续。

(2)将新的空间对象记录在叶子结点中：若L中尚有空位(即单元个数小于 M，M为一个结点中记录数目的最大值)，则在L中记录新的空间对象的 MBR(Minimum bounding rectangle最小边界矩形)和结点编号ID；若L中没有 空位，则将L结点中的单元(包括新的空间对象)分裂成两部分，分别记录在结 点L和新生成的结点LL中。

如何将一个裂隙圆盘构成的矩形集分裂成合适的两部分，是影响裂隙圆盘 检索效率的一个重要因素。本发明采用以面积作为标准算法将一个裂隙圆盘构 成的矩形集分裂成合适的两部分。即分裂后两部分的MBR的和最小。

(3)调整树的结构：从叶子结点L开始，对空间裂隙圆盘索引进行调整。

(3.1)初始化：令N为叶子结点L，NN为因新的空间对象插入而分裂生成 的新的结点；

(3.2)判断N是否为根结点：若N是根结点，则停止返回；

(3.3)调整父结点相应单元的I：调整结点N的父结点P中相应单元的I，使 其正好包含N结点中的所有单元；

(3.4)根据需要进一步分裂父结点：若存在分裂生成的新的结点NN，则在N 的父结点P中加入指向该结点的单元，若P结点中的单元的个数超过了M，则 需要分裂结点P，生成一个新的结点PP，并令N＝P，NN＝PP；转到步骤(3.2) 继续。

通过上面的分裂算法可以动态的加入空间对象，且不破坏空间裂隙圆盘索 引的平衡。

(4)生成新的根结点：若子结点的分裂导致根结点中单元个数超过M，则将 根结点分裂成两个结点，并再生成一个新的结点作为空间裂隙圆盘索引的根结 点，原根结点分裂成的两个结点作为其子结点；此时空间裂隙圆盘形成的索引 增加一层。

2.空间对象的删除操作：

(1)确定删除空间对象所在的叶子结点：通过空间裂隙圆盘索引的搜索算法， 确定删除空间对象所在的叶子结点L。

(2)删除空间对象：将结点L中的空间对象所在的单元删除。

(3)调整空间裂隙圆盘索引：从叶子结点L开始调整空间裂隙圆盘索引的操 作。这一步保证结点(根结点除外)中的单元的个数不低于其下限m。

(4)调整根结点：如果空间裂隙圆盘索引调整后，根结点只有一个孩子，则将根 结点删除，让其孩子结点成为空间裂隙圆盘索引的根结点，此时空间裂隙圆盘索 引减小一层。

由于空间对象的删除会导致结点单元的个数小于m，为了避免这种情况，提 高结点的利用率，需要对空间裂隙圆盘索引进行调整，保证每个结点(根结点除 外)的单元数不小于m。

调整算法包括以下步骤：

①初始化：令N为被删除单元的叶子结点L，Q为结点的集合，用以存放被 删除的结点，初始化为空；

②寻找N在父结点中的相应的单元：如果N为根结点就转向步骤⑥；反之 查找其父节点P，找到与结点N相对应的单元E_N；

③删除单元数小于m的结点：若结点N的单元数小于m，则在结点P中 删除单元E_N，并将结点N加入集合Q中；

④调整MBR：若结点N没有被删除，则调整单元E_N中的I，使其正好包 含结点N中的所有的单元；

⑤调整上一层：令N＝P，转到步骤②继续；

⑥被删除的结点重新插入到空间裂隙圆盘索引中：采用空间裂隙圆盘索引 的插入操作，将Q中的所有的被删除的结点重新的插入到空间裂隙圆盘索引中。

3.搜索算法：

空间裂隙圆盘索引的搜索算法是一个递归的过程。设搜索的区域为S，则搜 索区域S内的空间裂隙的过程如下：

(1)子树的搜索：从空间裂隙圆盘索引的根结点T开始，如果T不是叶子结 点，那么依次判断该结点中各单元的I与区域S的空间位置关系，若I与搜索 区域S相交，则该单元所指向的结点为待搜索的根结点，转到本步骤开始处继 续子树的搜索；若T是叶子结点则转至步骤(2)。

(2)叶子结点的搜索：若T是叶子结点，则依次判断其中的空间对象与搜索 区域S之间的空间位置关系，如果空间位置对象落在搜索区域S以内，那么表 明其满足搜索条件。

至此，完成了随机分布三维裂隙网络中裂隙对的初步相交判断，还可以通 过步骤五或者步骤六完成进一步的精确相交判断。

步骤五：将每块所有可能相交的裂隙对进行对接合并，去除其中重复的裂 隙对，对剩余可能相交的裂隙对进行精确相交判断，得到整个随机分布三维裂 隙网络中全部相交的裂隙对；其中精确判断具体为：根据已知裂隙圆盘的几何 参数，确定裂隙圆盘所在的平面，判断空间两个圆盘是否相交。

步骤六：对每块内所有可能相交的裂隙对进行精确相交判断，得到每块内 所有相交的裂隙对，将各块所有相交的裂隙对进行对接合并，去除其中重复的 裂隙对，得到整个随机分布三维裂隙网络中全部相交的裂隙对；其中精确判断 具体为：根据已知裂隙圆盘的几何参数，确定裂隙圆盘所在的平面，判断空间 两个圆盘是否相交。

精确判断时，应用步骤四中的方法，根据已知裂隙圆盘的几何参数，确定 裂隙圆盘所在的平面，判断空间两个圆盘是否相交。

以下介绍如何判断空间两个圆盘1和圆盘2是否相交。

根据已知裂隙圆盘1的几何参数(圆心坐标X₁、Y₁、Z₁，倾角α₁，倾向，β₁， 半径r₁)和圆盘2的几何参数(圆心坐标X₂、Y₂、Z₂，倾角α₂，倾向β₂，半径 r₂)，可知裂隙圆盘所在平面的方程为：

aX+bY+cZ+d＝0

其中，a＝Sinα₁×Sinβ₁，b＝Sinα₁×Cosβ₁，

c＝Cosα₁，d＝-(aX₁+bY₁+cZ₁)

两个裂隙圆盘所在的平面公式可表示为

$\left(\begin{array}{c}{a}_{1}X+b_{1}Y+c_{1}Z+d_1=0\\ a_{2}X+b_{2}Y+c_{2}Z+d_2=0\end{array}\right)$

两平面的夹角θ可表示为：

$\mathrm{Cos\theta }=\frac{a_1{a}_2+b_1{b}_2+c_1{c}_2}{\sqrt{(a_1^2+b_1^2+c_1^2)(a_2^2+b_2^2+c_2^2)}}$

θ≠0时即可求出两圆盘所在的平面的交线。

第i个裂隙圆心点到交线的距离为式：

$D_i=\sqrt{\frac{{((y_i-y_0)p-(z_i-z_0)n)}^2+{((z_i-z_0)m-(x_i-x_0)p)}^2+{((x_i-x_0)n-(y_i-y_0)m)}^2}{m^2+n^2+p^2}}$

$\left(\begin{array}{c}m=b_1{c}_2-b_2{c}_1\\ n=a_2{c}_1-a_1{c}_2\\ p=a_1{b}_2-a_2{b}_1\end{array}\right)$

其中x₀，y₀，z₀是交线上三维坐标中的任意点，m，n，p是交线的方向向量。

如果D₁+D₂＞r₁+r₂则不相交。否则，如果D₁+D₂≤r₁+r₂则需要通过圆心点 到两个平面交线的距离D₁、D₂与半径r₁、r₂的大小分别进行判断，只有当D₁＜r₁并且D₂＜r₂同时满足，两个圆盘才相交。

实施例：

实施本发明方法所用的计算机配置为Intel Core2 Quad CPU Q9500 (2.83＊2GHZ)，内存为4G，32位操作系统。不同的计算机配置将会对算法 的执行效率有所影响，本实施过程所采用的计算机配置为目前较为普遍机型。

根据本发明的随机分布三维裂隙网络中裂隙对相交判断方法包括以下步 骤：

(1)将空间中的每个裂隙用圆盘模拟。

将本发明的方法应用到9套随机分布三维裂隙网络系统中，各套裂隙系统 的几何参数见表1。选定的研究区域尺寸为50米＊50米＊50米，各套裂隙系统 几何参数体现出密度、直径(包括均值和方差)、产状(包括倾角和倾向的均 值和方差)以及组数的变化。其中第1，2，3，5，9套裂隙密度依次为0.5条 /立方米、0.6条/立方米、0.7条/立方米、0.8条/立方米、1.6条/立方米，以期 对比查看密度对算法的影响；第7套系统由2组裂隙组成，其余均为3组；第 4套裂隙系统的几何参数与第1套相比，其裂隙产状不同；第6套裂隙系统的 直径为3米，以期与第1套裂隙直径均值为2米做比较；第8套裂隙系统是与 第1套裂隙系统相比较，其直径的方差为2，小于第1套裂隙系统。

表1随机分布三维裂隙网络几何参数表

每个裂隙圆盘的(x，y，z，r，a，b)对应着几何参数的0≤x≤50米；0≤y≤50 米；0≤z≤50米；r取符合表中大小的均值和方差的随机数；a、b为表中产 状的倾角和倾向的随机数。

应用给定的三维随机分布裂隙几何参数，采用蒙特卡罗方法产生随机分布 的裂隙网络，同时产生三维空间中每个裂隙圆盘具体的三维几何参数；由于数 据量较大，以下显示第一套裂隙系统的前10条数据(共187500条裂隙数据)：

 50mx50mx50m   187500   x   y   z   r   a   b   1   17.2271   6.8811   4.1571   1.4513   18.4503   4.6561   2   7.2575   7.0847   9.0456   1.148   192.6735   2.5734   3   4.0107   18.1621   10.912   1.8614   269.2564   2.5406   4   1.9677   18.6152   4.4669   4.0442   163.2318   .4941   5   17.1971   9.9498   .0902   .1839   292.158   4.4939   6   2.9579   16.6138   3.1066   1.2771   264.507   3.8912   7   15.0786   18.3688   8.1482   12.3402   355.416   1.5901   8   19.2571   16.6133   19.5017   5.8991   329.7474   .762   9   18.8493   5.1012   19.7564   3.5206   269.2372   4.1167   10   2.4397   8.7232   14.7967   .1253   99.6888   4.8031

(2)用轴向包围盒代替空间中的裂隙圆盘以达到简化问题的目的；将10套 系统中全部的裂隙用轴向包围盒代替。分别计算三维空间中每个裂隙圆盘的外 接球体的外接正立方体的最小和最大坐标。由于数据量较大，以下显示第一套 裂隙系统的前10条裂隙的轴向包围盒最小和最大坐标数据(共187500条裂隙 数据)：

(3)将整个裂隙网络所在的岩块在空间中沿着三个坐标轴方向进行分块以 进一步降低问题规模；图3为两种分块方法示意图。

将每一个块作为独立单元，块内包含的裂隙和在切面上相交的裂隙归入该 块。

(4)块内采用空间数据索引算法，构建空间裂隙索引以快速进行相交判断得 到可能相交的裂隙对；以下显示第一套系统中与第1条裂隙相交的部分裂隙以 及相交的交线的终点坐标。

◆如下部分提供几种典型算法的优劣比较。

◆将裂隙网络在三维空间分块与不分块算法比较

理论上说，穷举在数据量小时有算法较为简单的优势，不需要花费额外的 时间和空间建立索引；而数据量较大时索引将大大减少访问时间。为了证明分 块是否能够大大降低问题的复杂度、块内穷举和索引有多大区别以及何种算法 最优。首先应用第1套裂隙系统的数据，分别采用不同分块方法及不同块内处 理算法进行测试，统计并对比运行时间。图5的第1和7套裂隙系统的数据在 分块与不分块算法情况下的运行时间，能够证明三种运行方案最佳时间中，不 分块全部穷举和索引方案的计算时间较长，采用了分块和索引的方法后对运算 时间有所改善。图中，

方案A：全部穷举算法

方案B：全部索引算法

方案C：分块后块内采用穷举和索引算法，根据不同的分块方式分为以下 四种具体方案：

C1：单方向分块+块内索引算法

C2：单方向分块+块内穷举算法

C3：三方向分块+块内索引算法

C4：三方向分块+块内穷举算法

◆基于分块的不同算法对比

图6是10套裂隙系统分别进行了三方向分块算法与单方向分块算法的模拟 结果。从图6可以看出：

①三方向分块方案比单方向分块方案效率高。

②相同分块大小情况下，单方向分块+穷举算法运算时间总是最多。

③在分块数目较大情况下，虽然三方向分块+穷举算法在分块数目较大时处 于整体算法比较的最低点，但是当分块数目不断减小，其运行时间也有增大趋 势，这意味着采用三方向分块+穷举算法时，分块的数目也有个限度，不是越多 越好。④同一套裂隙系统，同样分块大小，三方向分块+索引和三方向分块+穷 举算法的交叉点，处于块内圆心数目相对较少区域，而且三方向分块+穷举算法 占优势的区间相对较窄。

⑤三方向分块+穷举算法和单方向分块+索引算法对于分块的大小非常敏 感，尤其是三方向分块+穷举算法，当分块数目减少，块内裂隙数目增多时，分 块+穷举算法呈线性增长。而三方向分块+索引、单方向分块+穷举算法对分块 大小不太敏感。

也就是说，相比较单方向分块而言，三方向分块具有较好优势，而三方向 分块+块内索引算法具有稳定、可行性好等优势，是本发明所述有效算法。

(5)将每块所有可能相交的裂隙对进行对接合并，去除其中重复的裂隙对， 对剩余可能相交的裂隙对进行精确相交判断，得到整个随机分布三维裂隙网络 中全部相交的裂隙对。

以下显示第一套系统中与第1条裂隙相交的部分裂隙以及相交的交线的终 点坐标。

 裂隙1   与裂隙1相交的裂隙面序号   交线中点x   交线中点y   交线中点z   1   252   2.8119   16.5744   7.075   1   264   6.6466   15.7728   7.31   1   611   7.9653   17.8415   6.6089   1   766   3.5105   15.8494   7.3109   1   821   9.5689   17.276   6.7839

需要说明的是，可以先去除重复的裂隙对，再进行裂隙精确相交判断操作； 也可以先进行裂隙精确相交判断操作，再去除重复的裂隙对。本实施例采用的 是先去除重复的裂隙对，再进行裂隙精确相交判断操作。

另外，本发明中三维裂隙网络中裂隙对的初步或进一步精确相交判断后经 过分析处理将去除部分或全部对渗流路径无贡献的裂隙，从而有效提高渗流与 污染物运移数值计算效率。因为随机分布裂隙网络中存在长短不一裂隙，可能 包含较多对渗流无贡献裂隙，即使不进行进一步的精确判断，其初步判断结果 也可直接应用于识别渗流路径，其中可能含有一定数量未去除的对渗流无贡献 的裂隙，其对减少数值计算负荷与提高运算效率仍会具有实用意义。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知技术。

本发明不局限于权利要求和上述实施例所述及的内容，只要是根据本发明 的构思所创作出来的任何发明，都应归属于本发明的保护范围之内。