一种块体理论平衡区域图的三维可视化方法与流程

本发明涉及一种非连续工程计算可视化方法，特别是关于一种在岩土工程领域中应用的块体理论平衡区域图的三维可视化方法。

背景技术：

岩体作为边坡、地基和地下洞室等工程的直接作用对象，是在漫长的地质历史发展过程中经受了各种地质作用，并在地应力的长期作用下，在其内部保留了各种永久变形的迹象和地质构造形迹的天然地质体。由于岩体中存在各种结构面，结构面控制着岩体的稳定性，因此工程岩体与连续变形体相差甚远，一种计算分析方法能否在岩体稳定性分析中得到广泛应用，很大程度上取决于能否对岩体中的结构面进行合理模拟和计算。

作为岩体工程稳定性分析的一种重要方法，块体理论自石根华博士提出和建立以来，许多学者又相继从结构面的几何分布特征、物理特性以及力学行为等方面出发，并结合现代数学方法，判定岩石的稳定状况，正确认识受力岩体的变形和破坏规律，理论方面已经非常完备，并且在众多大型工程项目中得到应用。

倾角和倾向是表示空间平面几何方位的重要数据。倾角是岩层层面上的真倾斜线与其在水平面上投影线的夹角，表示在垂直断层面走向的直立剖面上该层面与水平面间的夹角；倾向是沿着岩层面倾斜方向向下引出垂直走向线的直线在水平面的投影地理方位，在空间直角坐标系中，x轴指向正东，y轴指向正北，z轴竖直向上。

结构面和临空面由倾角、倾向构建起来。岩体结构面分布于岩体的内部，块体理论把岩体结构面视为无限大平面，把岩体切割为不连续块体，由于结构面是块体失稳滑移的摩阻面，因此还需要考虑结构面的摩擦角以计算块体滑移时的摩擦阻力。临空面又分自然临空面和人工开挖面。临空面的存在是块体滑移的必要条件，因为块体在荷载的作用下只可能向临空面的方向即解除约束的方向移动。临空面和结构面一样，继承了抽象平面类的几何属性，但和结构面不同的是，临空面的几何方位在所研究的区域和阶段内是可变因素，或者说受外界条件影响较大，如不同的岩体工程结构就具有不同的临空面，尤其在工程前期选型阶段，通常要考虑多种不同的工程结构，也就是不同的临空面组合情况，通过分析计算选出最优的工程结构，而岩体结构面的几何方位或物理力学属性一般都是天然形成的，可视为常量。

块体是被各类结构面和临空面切割的岩体，亦称结构体。锥体是块体的数学抽象，它反映了块体表面的拓扑关系，是块体重要的属性之一；块体则是现实空间中存在的真正实体，有表面、边棱和顶点等几何要素。从块体理论的观点来看，非连续的块体系统是由块体组成，即块体一定是存在的，而其中的块体又可分为不同的类别。由结构面和临空面所切割的有限块体，如果它能够沿着不与其相邻块体相撞的一个方向滑动，那么这个有限块体是可动的。由结构面和临空面所切割的有限块体，如果它沿着任何方向滑动都会与其相邻块体相撞，那么这个有限块体是不可动的。

岩体的破坏可分为变形破坏和块体失稳破坏。岩体稳定性分析的目的就是通过各种手段和途径，正确认识受力岩体的变形和破坏规律，判定岩体的稳定状况，预测其未来的变化，制定出有效的工程处理措施。其一般研究过程包括：在地质勘查、岩体物理力学性质试验的基础上，通过经验类比、工程地质定性分析等，对岩体的稳定性及失稳破坏模式作定性分析；通过将复杂的岩体进行必要的简化，如抽象为刚性、弹塑性连续变形体，采用极限平衡分析法、有限差分法、有限单元法、离散单元法、边界单元法、非连续变形分析法和数值流行法等。

全空间赤平投影作为经典块体理论岩体稳定性分析方法的重要组成部分，其发展和应用最为成熟，在教学、科研、工程中发挥重要的作用。然而，由于赤平投影研究的是主动力合力不变的情况，对于如坝基和坝肩等复杂岩体稳定性的分析则显得尤为困难。为了解决上述问题，基于全空间赤平投影，通过选取单个节理锥结合平衡区域图以及滑动摩擦角等值线方法完美解决上述问题，正确认识受力岩体的变形和破坏规律，为工程规划、设计、施工和加固等工作提供了科学合理的建议和依据。虽然理论完备，但是平衡区域图发展至今，纵观国内外，尚未出现针对平衡区域图三维可视化方法分析的工具以及工程应用案例。

技术实现要素：

针对上述问题，本发明的目的是提供一种块体理论平衡区域图的三维可视化方法，该方法可以使分析结果更直观，更易于理解，能以多种方法演示同一结果，可使结果更具说服力。

为实现上述目的，本发明采取以下技术方案：一种块体理论平衡区域图的三维可视化方法，其特征在于包括以下步骤：1)建立二维平面投影区域以及与之对应的三维空间节理锥拾取对象模型；2)判断是否点击选择节理锥对象，选择则进入下一步，反之结束；3)求出节理锥各个平面指向外的法线方向，记录法线和相交棱信息；4)根据指向外的法线和相交棱信息，绘制各种运动模式下的平衡区域；5)计算小圆锥面起始向量和旋转角，绘制单面滑动摩擦角等值面；6)计算扇形平面的边界向量，绘制双面滑动摩擦角等值面；7)根据绘制的平衡区域图进行鼠标在三维空间的交互操作；8)计算主动合力所在的区域，判断出块体的运动模式：根据步骤5)和步骤6)分析块体达到稳定时，所需要达到的摩擦角值；当块体的实际摩擦角大于主动合力方向所在摩擦角等值面所对应的摩擦角值，块体保持稳定；反之，则为关键块体，需要提供支护措施。

进一步，所述步骤1)中，初始化节理锥对象的名字堆栈，绘制指定颜色下的全空间赤平投影图二维和三维选择拾取对象，并设置节理锥对应的名字。

进一步，所述步骤3)中，具体过程为：3.1)根据构成节理锥的每个结构面pi的二进制编码ai和法线向量计算出指向节理锥外的法线向量i＝1,……，n，n表示输入结构面的组数，第i组结构面pi由平移到球心的平面表示；3.2)依逆时针次序记录节理锥的各个平面的相交棱向量以及各个平面指向外的法线向量其中n表示输入结构面的组数，为平面和平面的相交棱；3.3)相交棱向量和构成平面为节理锥中平面指向外的法线向量。

进一步，所述步骤4)中，单面滑动平衡区域是由和以及相邻向量确定的平面围成的锥形区域，1≤i,j,k≤n，n表示输入结构面的组数；双面滑动平衡区域则是由3个向量和以及相邻向量确定的平面围成的锥形区域，其中1≤i≤n；稳定平衡区域是由指向节理锥外的法线向量以及由相邻法线向量确定的平面围成的锥形区域；掉落平衡区域则是由相交棱向量以及由相邻相交棱向量确定的平面围成的锥形区域。

进一步，绘制稳定和掉落平衡区域与参考球面相交的部分球面时，如果n＞3，则按照标号i1i2i3,i1i3i4,…,i1in-1in将向量围成的平衡区域划分为n-2个绘制区域；在使用3个向量绘制曲面时，计算出相邻向量间的中点并规范化至球面上，然后连接所有中点将原来的三角形划分为4个小三角形，记录各个小三角形的顶点向量；若对3个空间顶点向量进行m次迭代细分便能分解出4^m-1个小三角形，然后绘制拼接小三角形，逼近球面。

进一步，所述步骤5)中，单面滑动摩擦角等值面在三维空间中是圆锥点为球心o的小圆锥侧面一部分，该小圆锥侧面的绕轴为倾角为摩擦角φi，计算出与小圆锥侧面对应的单面滑动平衡区域的边界平面和的相交向量，该相交向量分别为小圆锥侧面的母线起始向量和母线终止向量母线起始向量和母线终止向量之间小圆锥侧面即为单面滑动摩擦角等值面。

进一步，所述单面滑动摩擦角等值面的绘制方法如下：5.1)计算母线起始向量绕旋转到母线终止向量时的旋转角γ；5.2)根据旋转角γ设置细分精度值为dγ，以为起始向量，递增旋转角，再根据罗德里格旋转矩阵计算出在小圆锥侧面的向量并保存到母线向量列表中，直到旋转夹角的值为γ；其中，还需要将加入到母线向量列表中；5.3)再用球心o与保存下的母线向量列表按照保存顺序的相邻两个向量组成小三角形，拼接成完整的小圆锥侧面，一个小圆锥侧面对应一个摩擦角值；5.4)根据摩擦角间距计算摩擦面的个数，绘制出该区域的所有单面滑动摩擦角等值面。

进一步，所述步骤6)中，双面滑动摩擦角等值面是以球心o为顶点，与双面滑动平衡区域的边界平面和分别交于母线起始向量和母线终止向量之间的扇形平面。

进一步，所述双面滑动摩擦角等值面的绘制方法为：6.1)将作为起始顶点向量并保存到顶点向量列表中；6.2)然后设置扇形平面细分迭代次数，每次迭代取中点并规范化，再保存到顶点向量列表中，最后将结束顶点向量加入到顶点向量列表中；6.3)用球心o与保存下的顶点向量列表按照保存顺序的相邻两个向量组成小三角形，拼接成完整的扇形平面。

本发明由于采取以上技术方案，其具有以下优点：1、本发明利用计算机图形学技术对经典块体理论的平衡区域图分析方法进行三维可视化，使得原本晦涩难懂的分析过程变得易于理解，分析结果更形象直观，从而促进关键块体理论在实际工程中的应用。2、本发明将倾角和倾向作为基本参数，将结构面视为无限大的平面，将块体视为凸体，将各种作用荷载视为空间向量，建立块体理论基本结构类型体系，进而研究块体的稳定性，进行岩体稳定性计算与工程开挖、支护分析，为工程设计提出科学合理的意见和建议。3、由于在实际工程中，特别是施工过程中，对岩石块体的稳定性分析通常采用关键块体理论，因此采用本发明的三维可视化方法展示其分析结果，一方面可以使分析结果更直观，更易于理解；另一方面结合全空间赤平投影、平衡区域、摩擦角等值线，以多种方法演示同一结果，可使结果更具说服力；同时，通过提供三维交互功能，可以对三维可视化结果进行某种交互式的测试，以验证分析结果的可靠性。

附图说明

图1是本发明的整体流程示意图；

图2是本发明的全空间赤平投影图二维和三维拾取对象示意图；

图3是本发明的单面滑动摩擦角等值面的生成示意图；

图4是本发明的双面滑动摩擦角等值面的生成示意图；

图5是本发明实施例中标注为0110节理锥的二维平衡区域图分析结果示意图；

图6是本发明实施例中标注为0110节理锥的各个运动模式下的三维平衡区域分析结果示意图；

图7是本发明实施例中标注为0110节理锥的三维单面滑动摩擦角等值面分析结果示意图；

图8是本发明实施例中标注为0110节理锥的三维双面滑动摩擦角等值面分析结果示意图；

图9是本发明实施例中标注为0110节理锥的完整的平衡区域图分析方法结果示意图。

具体实施方式

本发明基于块体理论三维球分析方法(jianyongli,jianxue,junxiao,yingwang.threedimensionalsphereanalysismethodofblocktheory,2010internationalconferenceoncomputerapplicationandsystemmodeling(iccasm2010),taiyuan,shanxi,pages578-582,2010.10.)建立的节理锥(jp)模型，在节理锥模型的基础上提出一种块体理论平衡区域图的三维可视化方法，下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明提供一种块体理论平衡区域图的三维可视化方法，将分析过程和结果映射到三维空间，同时满足交互的要求，其关键是根据选定的节理锥(jp)，用三维图形的方式直接表示其平衡区域和摩擦角等值面，并给出滑动模式以助于判断块体稳定性。该方法所需的输入信息包括：n组结构面的倾角(dip)、倾向(dipdirection)和摩擦角(frictionangle)；m组临空面的倾角(dip)和倾向(dipdirection)；其包括以下步骤：

1)如图2所示，建立二维平面投影区域以及与之对应的三维空间节理锥拾取对象模型：初始化节理锥对象的名字堆栈，绘制指定颜色下的全空间赤平投影图二维和三维选择拾取对象，并设置节理锥对应的名字。

其中，采用opengl(opengraphicslibrary，开放图形库)进行绘制。

2)判断是否点击选择节理锥对象，选择则进入下一步，反之结束；

当点击全空间赤平投影的各个不同颜色区域时，opengl反馈对应的节理锥名字，然后返回节理锥对象。

3)求出节理锥各个平面指向外的法线方向，记录法线和相交棱信息：

3.1)根据构成jp的每个结构面pi的二进制编码ai和法线向量计算出指向节理锥外的法线向量可用下式计算：

式中，i＝1,……，n，n表示输入结构面的组数，第i组结构面pi由平移到球心的平面表示。每个结构面将整个空间分为两个半空间，jp实际上是由各结构面所取的半空间相交而成的，按每一结构面上下2个半空间的组合，n个结构面最多可能围成2ⁿ个jp，二进制编码ai是用来指定取第i个结构面的哪个半空间，0表示上半空间(即法向量指向的半空间)，1表示下半空间，则所有结构面上下半空间的组合可表示为从0到2ⁿ-1的二进制数，每个二进制数对应一个jp。

3.2)依逆时针次序记录节理锥的各个平面的相交棱向量以及各个平面指向外的法线向量为平面和平面的相交棱；k＝1,……，n。

3.3)相交棱向量和可以构成平面为jp中平面指向外的法线向量。

4)根据指向外的法线和相交棱信息，绘制各种运动模式下的平衡区域：根据法线向量、相交棱向量之间的关系，结合opengl图形库可绘制出稳定、掉落、单面滑动、双面滑动模式下的图形区域；

单面滑动平衡区域是由3个向量和以及相邻向量确定的平面围成的锥形区域，其中1≤i,j,k≤n，n组结构面可形成n个单面滑动平衡区域，即单面滑动平衡区域个数和结构面组数相等。

同理，双面滑动平衡区域则是由3个向量和以及相邻向量确定的平面围成的锥形区域，其中1≤i,j≤n，n组结构面可形成n个双面滑动平衡区域，即双面滑动平衡区域个数和结构面组数相等。

稳定平衡区域是由指向节理锥外的法线向量以及由相邻法线向量确定的平面围成的锥形区域。

掉落平衡区域则是由相交棱向量以及由相邻相交棱向量确定的平面围成的锥形区域。

在绘制稳定和掉落平衡区域与参考球面相交的部分球面时，如果n＞3，则可以按照标号i1i2i3,i1i3i4,…,i1in-1in将向量围成的平衡区域划分为n-2个绘制区域。需要注意的是，在使用3个向量绘制曲面时，计算出相邻向量间的中点并规范化至球面上，然后连接所有中点将原来的三角形划分为4个小三角形，记录各个小三角形的顶点向量。若对3个空间顶点向量进行m次迭代细分便可分解出4^m-1个小三角形，然后绘制拼接小三角形，m越大则绘制的图形越逼近球面。

5)计算小圆锥面起始向量和旋转角，绘制单面滑动摩擦角等值面：

单面滑动摩擦角等值面在三维空间中是圆锥点为球心o的小圆锥侧面一部分，该小圆锥侧面的绕轴为倾角为摩擦角φi，计算出与小圆锥侧面对应的单面滑动平衡区域的边界平面和的相交向量(边界平面表示由向量和向量确定的平面，边界平面含义类似)，该相交向量分别为小圆锥侧面的母线起始向量和母线终止向量母线起始向量和母线终止向量之间小圆锥侧面即为单面滑动摩擦角等值面；其中：

式中，摩擦角φi的取值范围为0度到90度，每个给定φi对应一个摩擦角等值面，即φi从0到90度取不同的值可绘制出一组不同的等值面。

5.1)计算母线起始向量绕旋转到母线终止向量时的旋转角γ；其中，母线从位置绕旋转到位置扫过的区域即为所需绘制的曲面；

5.2)根据旋转角γ设置细分精度值为dγ，以为起始向量，递增旋转角，再根据罗德里格旋转矩阵计算出在小圆锥侧面的向量并保存到母线向量列表中，直到旋转夹角的值为γ；其中，由于精度值dγ的计算误差可能会导致旋转角γ后得到的向量与不相等，故还需要将加入到母线向量列表中；

5.3)再用球心o与保存下的母线向量列表按照保存顺序的相邻两个向量组成小三角形，拼接成完整的小圆锥侧面，一个小圆锥侧面对应一个摩擦角值；

其中，母线向量起点(球心o)和两条母线向量各自的终点，由这3点可构成一个三角形，每两条相邻的母线向量即可确定一个三角形，由这些相邻三角形最终拼接成小圆锥侧面，如图3所示。

5.4)根据摩擦角间距计算摩擦面的个数，绘制出该区域的所有单面滑动摩擦角等值面。

需要注意的是单面滑动摩擦角等值面与单面滑动平衡区域相对应，可在各个单面滑动平衡区域上显示单面滑动摩擦角等值面的变化情况。

6)计算扇形平面的边界向量，绘制双面滑动摩擦角等值面：双面滑动摩擦角等值面是以球心o为顶点，与双面滑动平衡区域的边界平面和分别交于和之间的扇形平面，其中和的计算公式如下：

6.1)将作为起始顶点向量并保存到顶点向量列表中；

6.2)然后设置扇形平面细分迭代次数，每次迭代取中点并规范化，再保存到顶点向量列表中，最后将结束顶点向量加入到顶点向量列表中。

6.3)用球心o与保存下的顶点向量列表按照保存顺序的相邻两个向量组成小三角形，拼接成完整的扇形平面(如图4所示)。

需要注意的是双面滑动摩擦角等值面与双面滑动平衡区域相对应，可在各个单面滑动平衡区域上显示单面滑动摩擦角等值面的变化情况。

7)平衡区域图三维可视化结果采用opengl图形库进行绘制，并根据平衡区域图三维可视化结果进行鼠标在三维空间的交互操作：

平衡区域图三维可视化结果采用opengl图形库进行绘制，以适应跨平台的需求，而三维交互操作借助鼠标或其他指点设备(例如触摸屏等)完成，同时要求交互的反馈与直观经验一致，这里涉及两个需要处理的问题，一个是对绘制结果中的三维图形对象的选取，这一点借助opengl选择模式下的绘制来完成；另一个是将鼠标或二维指点设备的运动转换为三维空间中被选中的图形对象的运动，由于鼠标的运动只能给出二维的坐标信息，无法直接恢复出三维坐标，需要其他辅助信息的加入才能完成转换，因此将鼠标的移动限制在屏幕坐标系的xoy平面(即显示屏幕所在平面)即可得到鼠标移动的第三个坐标分量z＝0，当转换得到鼠标移动在对象坐标系中三维向量后，还需将其投影到jp的滑动方向上，以使jp只能沿其滑动方向移动，从而使交互操作更具实际意义。上述二维运动到三维运动的转换可通过下式的计算完成：

v＝((m^-1v^-1p^-1u)·s)s(7)

其中，m为图形引擎的模型变换，v为图形引擎的视图变换，p为图形引擎的投影变换矩阵，u是鼠标在归一化设备坐标系(ndc)下的移动向量，s为jp的滑动方向向量；

8)计算结果分析：计算主动合力所在的区域，判断出块体的运动模式。

根据步骤5)和步骤6)分析块体达到稳定时，所需要达到的摩擦角值；当块体的实际摩擦角大于主动合力方向所在摩擦角等值面所对应的摩擦角值，块体保持稳定；反之，则为关键块体，需要提供支护措施。

综上所述，本发明通过上述步骤，即可实现平衡区域图分析的三维可视化，为验证本发明所提方法的有效性和实用性，下面给出一个实例。表格1为验证实例的输入数据，其中，倾角、倾向确定平面朝向(法向量方向)，其二维平衡区域图分析结果如图5所示，通过采用本发明的方法所得的三维可视化结果如图6～图9所示，其中，图9为整体可视化的结果，图6为标号0110的jp稳定、掉落、单面滑动以及双面滑动平衡区域的可视化结果，图7为标号0110的单面滑动摩擦角等值面可视化结果，图8为标号0110的双面滑动摩擦角等值面可视化结果。

表1验证实例的结构面、临空面输入数据

上述各实施例仅用于说明本发明，各部件的结构、尺寸、设置位置及形状都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别部件进行的改进和等同变换，均不应排除在本发明的保护范围之外。