2.4 三维块体几何识别算法

2.4.1 几何识别算法的基本思路

由2.2节可知，任意形状的空间多面体可以看作三维空间中一有向复合形K，此复合形必须满足式（2.7）。由此可知，对于多面体的每一条边在空间上都与多面体两个环路的两条有向边相重合，而这两条有向边的起点和终点在空间上位置正好颠倒，因此其和为零。在此结论基础上，可得到三维块体几何识别算法的基本思路，即通过裂隙面之间求交获取点，进而得到边，经边的正则化（边的“树枝删除”）处理之后，利用边界算子并按相应规则搜索形成空间环路（即二维块体），再经面的正则化（面的“树枝删除”）处理之后利用边界算子并按相应规则进行搜索，最终形成三维块体。环路和块体的搜索是一个反复利用边界算子的过程，直至满足∂（∂K）=0才结束，相应的搜索规则可采用最大右转角准则^［33］或最小左转角准则^［147］。

图2.10给出了利用边界算子搜索形成平面环路（即二维块体）的过程^［136］。

图2.10 二维块体生成中的边界算子运算

（a）形成的二维块体；（b）顶点标记；（c）有向边的分解和端点对的形成；（d）～（i）搜索流程

其中，（d） 以边（a⁺，b^-）为起始边，∂（a，b）=a⁺+b^-；

（e） 增加边（b⁺，c^-），∂［（a，b）+（b，c）］=a⁺+b^-+b⁺+c^-=a⁺+c^-；

（f） 增加边（c⁺，d^-），∂［（a，b）+（b，c）+（c，d）］=a⁺+c^-+c⁺+d^-=a⁺+d^-；

（g） 增加边（d⁺，e^-），∂［（a，b）+（b，c）+（c，d）+（d，e）］=a⁺+d^-+d⁺+e^-=a⁺+e^-；

（h） 增加边（e⁺，f^-），∂［（a，b）+（b，c）+（c，d）+（d，e）+（e，f）］=a⁺+e^-+e⁺+f^-=a⁺+f^-；

（i） 增加边（f⁺，a^-），∂［（a，b）+（b，c）+（c，d）+（d，e）+（e，f）+（f，a）］=0，即完成二维块体搜索。

图2.11 三维块体生成中的边界算子运算

在完成二维块体搜索之后，可在所形成的面环基础上利用边界算子进行三维块体的搜索^［137］。图2.11给出这一过程，即

（a）以面2为开始面，∂∂K=-（e₁）+（e₅）+（e₆）+（e₇）；

（b）增加面6，∂∂K=-（e₁）+（e₅）+（e₁₀）+（e₉）+（e₈）+（e₇）；

（c）增加面5，∂∂K=-（e₁）-（e₄）+（e₁₂）+（e₉）+（e₈）+（e₇）；

（d）增加面1，∂∂K=（e₂）+（e₃）+（e₁₂）+（e₉）+（e₈）+（e₇）；

（e）增加面3，∂∂K=（e₃）+（e₁₂）+（e₉）+（e₁₁）；

（f）增加面4，∂∂K=0，三维块体搜索结束。

根据以上分析，可得出三维块体几何识别算法具体实现流程，如图2.12所示。

2.4.2 裂隙面表示及迹线生成

对于三维块体的自动识别问题，首先需要求解构成空间块体棱边的裂隙面相交迹线。假定裂隙面是平面且光滑，具有有限尺寸的空间多边形；裂隙面可通过实测地质资料生成，也可以由随机过程产生。

1.裂隙面表示

为了表征三维空间中的裂隙面，除了裂隙面的倾向ω和倾角ρ以外，还需要裂隙平面内的参考点 （如质心或几何中心等）和描述裂隙面形态及尺寸的参数，以上这些参数中，裂隙面形态及尺寸参数的确定是个难点。目前，对裂隙形状主要有两种假设：圆盘形和多边形。Chile认为，在岩体完整性较好、裂隙发育较小的情况下，圆盘模型有一定的适应性；当岩体中大裂隙较多、岩体被切割成块状时，则多边形模型较合理；Dershowitz通过对野外裂隙的分析认为，多边形模型较符合实际情况。因此，在研究三维块体几何识别问题时，裂隙面形态采用多边形假定。

图2.12 三维块体几何识别流程图

为了分析和处理上比较方便，同时又能考虑其他裂隙面的截断作用，本书假定裂隙面的初始形态为内接于圆形的正多边形，若其延伸范围受到其他裂隙面的约束，则修正为一般多边形。因此，本书采用一个六参数d_i（x_i，y_i，z_i，ω_i，ρ_i，r_i）集合和一组有限数目且拥有具体空间坐标值的角点集合来联合表示空间的裂隙面，此角点集合为V（X）={v_i（x_i，y_i，z_i），i=1，2，…，N_v，（x_i，y_i，z_i）∈R³}，N_v表示裂隙面的角点总数，所有这些角点必须与裂隙平面共面。需要说明的是：代表裂隙面的空间多边形均为凸多边形。

因此，裂隙面多边形模型的生成包括两步：即正多边形的生成和一般多边形的修正。

为了便于生成内接于圆形的正多边形，假定裂隙面拥有两套坐标系，即整体坐标系O—XYZ和局部坐标系o-n-s-t，具体如图2.13所示。

其中，图2.13（a）定义了裂隙的局部坐标系o-n-s-t，满足右手系规则，n轴正指向裂隙面的上盘方向，s轴沿地面走向，t轴沿倾角方向；图2.13（b）反映了裂隙［x₁（t）-y₁（s）-z₁（n）］的整体坐标系与局部坐标系的转换。

图2.13 整体与局部坐标系

假设第i个裂隙面中心整体坐标为，此裂隙面的整体坐标系与局部坐标系转换见式 （2.15）。其中，为局部坐标系与整体坐标系的转换矩阵，此矩阵由正北方向与X轴夹角θ（与整体坐标系的选择有关，与裂隙面无关）、第i个裂隙面倾向ωⁱ和倾角ρⁱ决定；矩阵Tⁱ是正交单位矩阵，即

假定初始多边形为内接正十二多边形（图2.14），则此多边形的角点局部坐标可由式（2.17）求得，其整体坐标可通过式（2.18）转换而得。

图2.14 裂隙面的初始形态

在形成初始正多边形各个角点后，需要检查裂隙面的受约束情况，若受到其他裂隙面的截断作用，则需要对其进行修正，转化为几何图形学语言就是对此裂隙面进行面-面相减的布尔操作进而选取符合约束情况的裂隙面片部分；若没有受到约束，则直接将初始多边形转化为一般多变形。经过转化处理以后，便可形成裂隙面形状及尺寸的角点集合表示，即式（2.19）。

2.迹线生成

裂隙面之间或裂隙面与边界面之间的迹线生成问题，从几何图形学角度来看，就是一般多边形之间的求交问题。

令描述第i和j个裂隙面信息的参数集合和角点集合为：

其所在平面的法向矢量为），平面方程为：

若满足式（2.26），则裂隙面i和j平行或重合；否则，两者所在平面相交。

令，则两平面的交线方程为：

为便于计算，将描述第i个裂隙面形态及尺寸的角点集合［见式（2.21）］可转化为如下线段集合形式，见式（2.28）。

将式（2.27）与式（2.28）联立求解，就可得到两裂隙面所在平面相交直线与裂隙面i的交点。由于裂隙面i为凸多边形，所产生的交点不会超过两个；若存在两个交点即（x₁，y₁，z₁）和（x₂，y₂，z₂），则表明第i裂隙面与第j裂隙面所在的平面相交；若不存在，则表明这两个裂隙面不相交。

同样，对第j裂隙面进行上述运算，可得到第j裂隙面与第i裂隙面所在的平面的交点，即（x₃，y₃，z₃）和（x₄，y₄，z₄）；或者两平面不相交。

这样便产生了四个交点P_i（x_i，y_i，z_i），i=1，2，3，4，它们在整体坐标下是共线的，裂隙面相交产生的实际迹线由线段的公共相交部分决定 （图2.15）。

在完成迹线生成之后，需要引入一个矩阵用来记录裂隙面之间的连接（相交）关系。假设裂隙总数为N_d，那么表示这种连接关系的矩阵为一N_d×N_d方阵，即

矩阵C^d的元素为

图2.15 两条共线的相交迹线段之间求交

对于矩阵C^d的对角元素，有

矩阵C^d是一个对称矩阵，被称作裂隙面连接矩阵。矩阵C^d第i列所有元素之和，表示与第i裂隙面相交的所有裂隙面总数。

如果采用常规储存技术，那么矩阵C^d将是非常大的。因此，当考虑大量裂隙面的实际工程问题时，需要采用压缩储存技术或链表结构来记录和检索裂隙间的连接关系，这对渗流分析和块体识别来说是非常重要的。

3.裂隙面与迹线信息存储

在完成迹线生成之后，输入的和生成的数据集合有：

（1）裂隙面的棱边数据：D（E）=d{（e_i，1，e_i，2，…，e_i，M+K），i=1，2，…，N_d}，在所有棱边e_i，j（1，2，…，M+K）中，定义裂隙面几何形态的边界棱边数为M个，定义该裂隙面与其他裂隙面的相交迹线的内部棱边数为K个。

（2）角点坐标数据：V（X）={（x_i，y_i，z_i），i=1，2，…，N_v}，其中，N_v为当前角点总数，（x_i，y_i，z_i）是角点的整体坐标，这些角点既包括定义所有裂隙面几何边界的顶点，也包括所有的相交迹线的端点。

（3）棱边数据：E（V）={（v_i，1，v_i，2），i=1，2，…，N_e}，其中，N_e为当前棱边总数，既包括定义裂隙面边界的棱边，也包括裂隙面间的相交迹线。

（4）裂隙面的相交数据：当前裂隙面连接（相交）矩阵C^d或与之等价的链表数据结构。

2.4.3 空间面环的几何识别

空间环路识别部分主要包括：首先是求解每个裂隙面上所有棱边（即相交棱边和边界棱边）间的交点，进而删除每个裂隙面上的孤立和“悬空”的迹线段（即棱边的正则化过程），然后对每个裂隙进行空间环路搜索，判别环路之间隶属关系，最终生成裂隙面上的有向面。

1.迹线相交

假设是某个裂隙面上的第i和j个相交棱边的起始点和终止点，如图2.16 （a）所示。如果这两者相交，将产生新的交点V（x_s，y_s，z_s），如图2.16 （b）所示。

图2.16 两条迹线段相交

这两个迹线段的参数方程如下：

联立式（2.29）和式（2.30），可得

若令求解式（2.29），可得

如果同时满足0≤t_i≤1，0≤t_j≤1，则表明这两个棱边相交，将t_i代入式（2.29）或t_j代入式（2.30）便可产生新交点（x_s，y_s，z_s）；否则，这两条棱边不相交。对包括裂隙面边界棱边在内的所有的相交棱边执行上述运算，便可完成迹线相交。

在完成迹线求交之后，为了棱边正则化的方便，需要对迹线段所在直线上的所有交点（包括端点）进行排序，排序可采用多种算法，如冒泡算法、快速排序法等。

2.边的正则化（又称边的“树枝删除”）

以上形成的交点和棱边网络能够形成一个完全连通的、没有“树枝边缘”的、可被看做是平面复形的二维图结构。因此，对于每个裂隙面的交点和棱边需满足以下条件：①每个角点至少与两条棱边相连；②每条棱边由两个角点组成且是两个面的公共交线；③最少有三个角点和棱边。

在进行边的正则化时，那些“悬空”的迹线段被删除后，有些不在裂隙面边界上的端点将变成孤立角点，这些孤立点也要被删除，这是一个反复的过程，直至最后形成“无树枝”的有向图。“删除树枝”意味要改变裂隙面的相交关系，因此，需要随时修改裂隙面连接矩阵C^d中的相关元素。

与此同时，需随时更新棱边数据E（V）、角点数据V（X）以及裂隙面的棱边数据D（E）。此外还需对每个裂隙面建立两个新的数据集合，一个是棱边线段矩阵S（V_i）={（v_i，1，v_i，2，…，v_i，k）_j，i=1，2，…，N_s；j=1，2，…，N_d}_j，用于记录第j裂隙面上每个迹线段所在直线上所有交点（包括端点）的序列；另外一个就是角点-棱边连接矩阵V（E）_j={（e_i，1，e_i，2，…，e_i，k）_j，i=1，2，…，N_v；j=1，2，…，N_d}，用于记录拥有公共角点i的所有棱边。

另外，为了便于边的正则化运算，还需要建立每个裂隙面的棱边连接矩阵，是此裂隙面上沿某条迹线段上的最大交点数。矩阵的每列元素v_k1，v_k2，…，v_kn是根据距离相应迹线段起始端点的远近而确定的交点序列，此矩阵反映了此裂隙面上迹线段之间的相交关系。交点v_kl出现在矩阵的第k行和第l列，这表明它是迹线段k和l的公共点。

在完成“树枝删除”后，需要对形成的棱边进行定向操作。故引入辅助矩阵O=［o_ij］，用于标记面环搜索中棱边方向的使用状态。此矩阵称为棱边方向索引矩阵，矩阵元素如下：

o_ij=2 （棱边的初始状态，两个方向均没有使用）；

o_ij=1 （棱边的中间状态，从角点v_i，j到角点v_i，j+1的方向没有使用）；

o_ij=-1 （棱边的中间状态，从角点v_i，j+1到角点v_i，j的方向没有使用）；

o_ij=0 （棱边的最终状态，两个方向都被使用）。

图2.17给出了一个含有19条相交迹线和25个角点的裂隙面及相应的棱边线段矩阵。

图2.17 裂隙面上的相交迹线及相应棱边表示

（a）迹线段及顶点；（b）棱边线段矩阵S（V）

1～19—相交迹线编号；·—角点；①—角点编号

图2.18给出了图2.17中裂隙面在完成“树枝删除”后的情况。关于“树枝删除”的详细过程可参见石根华的著作^［147］。矩阵O的所有元素均为2，在完成空间环路识别后，矩阵中所有元素值为0。

图2.18 “树枝删除”后的裂隙面表示

（a）无树枝的有向图；（b）棱边线段矩阵S（V）；（c）棱边方向索引矩阵O

3.空间环路搜索及有向面的生成

空间环路搜索是从经正则化处理后的裂隙面上任一角点及与之相连的任一棱边开始，共用此角点的所有棱边均可从角点—棱边连接矩阵V（E）_j中获取。以当前棱边为起始边，围绕当前角点反时针旋转进而寻找与起始边构成最大右转角的棱边，并把它作为下一条棱边进行搜索，这就是最大右转角准则。由于棱边具有两个方向，因此，当每条棱边均被搜索两次时，环路搜索才算结束。

图2.19 环路搜索准则示意图

需要说明的是：最大右转角准则和最小左转角准则本质上相同的，只是转动方向不一样而已。

最大右转角准则是环路搜索的重要判别准则，如图2.19所示，设e₀边所在环路下一条可能的有向边为e_i（i=1，2，…，n），设e_i对应的右转角为θ_i：

式中：e^₀、e^_i为e₀、e_i对应的有向边的矢量；n^为有向边所在有向面的法向矢量。而θ_max=max（θ₁，θ₂，…，θ_n）=θ_k，则θ_max所对应的有向边e_k即为f的下一条环路有向边。

按照上述准则会产生两种类型环路：有内侧域的环路 （即内环）和有外侧域的环路（即外环），具体如图2.20所示。环路的类型可以根据环路的有向面积A的正负来判断，当A＞0时，为内环，否则为外环。环路有向面积的计算可参看2.6节一般多边形的面积计算式 （2.87）。

图2.20 环路拓扑类型

（a）内环；（b）外环

由于面是由一个且仅有一个外环和若干个内环组成的，因此，需要对所形成的空间环路进行包含隶属关系检查，以确定有向面的环路组成。定义空间环路的包含关系如下：若环1所包含有界平面，沿环2前进时其左手方向指向无限域，两环实际所围成平面有公共部分，则环1包含环2且环1为外环，环2为内环。由环搜索的最大右转角准则已经保证各环路不相交亦不自交，因此环的包含关系的检测就简化为点与面之间包含关系性的检测。

对于任意环路L：①若检测发现有环路包含其内，则环路L为某面的外环，与其余所有包含在此外环内的内环一起形成该面的环路；②若检测发现有其他反向环路包括环路L，则此环路为内环，相应外环L与内环组成该面的环路；③若检测发现环路L与其他环路无包含关系，则此环路为外环并为单独形成有向面的环路。图2.21给出了空间环路搜索及有向面生成的例子。若以棱边e_1，6为起始边，角点v₁为起始角点，角点v₆为特征角点，那么围绕角点v₆的边有棱边e_6，5和e_6，15；从图2.21（a）中可以看出，棱边e_1，6与e_6，5形成的右转角要小于棱边e_1，6与e_6，15之间的右转角，按照最大右转角准则，选取棱边e_6，15为检索棱边以及角点v₁₅为特征角点继续进行搜索；接着搜索到棱边e_15，4，特征角点为v₄；进而搜索到棱边e_4，1，此时特征角点为v₁，这与起始角点重合，说明已形成一个环路；而后需要判断环路类型，由于环路的点序列为v₁→v₆→v₁₅→v₄→v₁，呈逆时针旋转，所形成的有向面积大于0，故为内环1号。接着继续搜索形成2号环路，其点序列为v₁→v₄→v₃→v₂→v₁₀→v₉→v₁₁→v₁₂→v₇→v₈→v₅→v₆→v₁，呈顺时针旋转，所形成的有向面积小于0，故为外部环路。按照这种方法可以搜索到内部环路3号～10号以及一个与10号内环角点完全相同的外部环路，具体见图2.21（c）～（j）；接着需要判断这些环路之间相互关系，进而形成有向面环。这就形成了有向面环1号、3号～10号，其中8号有向面环包括两个环路［图2.21（k）］。

图2.21（一） 环路搜索示例

（a）形成二维块体1；（b）形成二维块体2—外部块体；（c）形成二维块体3；（d）形成二维块体4；（e）形成二维块体5；（f）形成二维块体6；（g）形成二维块体7；（h）形成二维块体8—包含二维块体10

图2.21（二） 环路搜索示例

（i）形成二维块体9；（j）形成二维块体10；（k）形成二维块体8的内环

2.4.4 空间块体的几何识别

空间块体识别部分主要包括：在空间面环搜索的基础上，删除那些孤立的或“悬空”的有向面（即面的正则化过程），然后按照一定准则进行空间块体搜索和判断，进而生成空间块体。

1.面的正则化（又称面的“树枝删除”）

由于三角形是2维单纯形，因此，面的正则化过程首先从移去或删除那些拥有与其他面共用的非共线边数目少于三条的面开始；如果某个面的一条边或多条边不能被其他面所共用，那么这个面就是孤立面，应当给予删除；当删除一些面后，其他一些面将失去它们的公共边而变成“悬空”面，这些面也应当删除。这是一个反复的过程，直至形成连通的、没有“树枝边缘”的空间裂隙面片网络。

与此同时，还需建立两个非常有用的矩阵，即棱边-环路连接矩阵和环路—环路连接矩阵，如下

式中：m为所有共用i边的有向面总数；n为与有向面i相交的有向面总数；N_e为总的棱边数；N_f为总的有向面数；f_i，j用于记录有向面的状态。

在完成“树枝删除”后，需要对形成的有向面进行定向操作。故引入辅助数组{I_F}={（I_i），i=1，2，…，N_F}，用于标记块体搜索中有向面方向的使用状态。此数组称为有向面方向索引矩阵，数组元素如下：

I_i=2 （有向面i的初始状态，其两个方向均没有使用）；

I_i=1 （有向面i的中间状态，其负方向被使用）；

I_i=-1 （有向面i的中间状态，其正方向被使用）；

I_i=0 （有向面i的最终状态，其两个方向均被使用）。

2.空间块体搜索及有向体的生成

有向面的搜索本质上是二维问题，而有向体的搜索是三维问题，该问题可以看成是对于二维问题的推广。有向边的方向为其实际空间走向，而有向面的方向则为其外法线矢量方向，体搜索的准则是二维环路搜索准则的推广。

图2.22 块体搜索准则示意图

空间块体搜索选择某一有向面作为起始环路，选择此有向面上某一条边为起始棱边，围绕当前边反时针旋转进而寻找与起始有向面构成最大右转角的环路，并把它作为下一个有向面，并选择此有向面其他棱边作为检索棱边继续进行搜索。由于有向面和棱边具有两个方向，因此，当所有的有向面和棱边均被检索两次，块体搜索才算结束。

图2.22表示体（壳）搜索的最大右转角准则。对于空间已搜索到的有向面f₀，确定以f₀的某边e为公共边的下一条有向边，设满足与f₀有公共边e的有向面f_i（i=1，2，…，n），f_i与f₀的右转角为：

式中：n^₀为f₀外法向矢量；n^_i为f_i外方向矢量。

则θ_max=max（θ₁，θ₂，…，θ_n）=θ_k所对应的有向面f_k（1≤k≤n）为相应的下一个有向面。

按照上述准则会产生两种类型的块体：有内侧域的块体（即内壳）和有外侧域的块体（即外壳），具体如图2.23所示。块体的类型可以根据块体的有向体积V的正负来判断，当V>0时，为内壳，否则为外壳。有向体积的计算参看2.6节一般多面体的体积计算式（2.93）。

由于有向体是由一个且仅有一个外壳和若干个内壳组成的。因此，需要对所形成的块体进行包含隶属关系检查，以确定有向体的块体组成。由体 （壳）搜索的最大右转角准则已经保证各体 （壳）不相交亦不自交，因此块体的包含关系的检测就简化为点与体之间包含关系性的检测。

图2.23 块体拓扑类型

（a）外壳；（b）内壳

图2.24～图2.29给出空间块体搜索例子，其实现流程与空间面环搜索过程相同。

图2.24 空间块体B₀的搜索

（a）块体B₀及相应识别信息；（b）由处于环路群中最外边界处环路所形成的外部块体

图2.25 空间块体B₁的搜索

（a）块体B₁及识别信息；（b）块体B₁单独表示

图2.26 空间块体B₂的搜索

（a）块体B₂及识别信息；（b）块体B₁和B₂单独表示

图2.27 空间块体B₃的搜索

（a）块体B₃及识别信息；（b）块体B₁、B₂和B₃单独表示

图2.28 空间块体B₄的搜索

（a）块体B₄及识别信息，带十个面具有复连通域性质的块体；（b）块体B₁、B₂、B₃和B₄单独表示；（c）未经块体搜索的面环

图2.29 空间块体B₅的搜索

（a）添加最后的环路到B₅；（b）块体B₅识别；（c）搜索形成的五个块体