CN103759677B - 基于三角面积比法度量岩体结构面三维粗糙度的方法 - Google Patents

Abstract

本发明一种基于三角面积比法度量岩体结构面三维粗糙度的方法，属于土木工程中岩土技术领域，本发明应用数字摄影测量系统进行结构面的采集，通过结构面出露部分的表面形态来反映整个结构面的三维空间形态，通过最小二乘法对结构面进行平面拟合，计算结构面三角网格的总面积与平面面积之比，反应了拟合结构面的平整程度；本发明充分利用了摄影测量法可以快速准确的测量结构面三维数据的优势，在此基础上应用Delaunay三角剖分方法进行结构面表面的重建，实现通过三维三角网格来反映结构面的表面形态；在结构面表面三维重建的基础上，提出了比表面积法，通过结构面的网格面积与拟合面积之间的比值来反映粗糙度大小。

Description

基于三角面积比法度量岩体结构面三维粗糙度的方法

技术领域

本发明属于土木工程中岩土技术领域，具体涉及一种基于三角面积比法度量岩体结构面三维粗糙度的方法。

背景技术

岩体是由岩块和分离切割岩块的结构面组成的岩体结构。岩体结构面具有大小不一、产状不同、形态各异等特点，因而岩体呈现出不连续、非均质及各向异性等力学性质。大量的工程实践可以表明在隧道洞室开挖过程中，结构面的几何和物理特性对岩体的变形和破坏起着关键性的控制作用，对工程岩体的稳定性有重要影响。

岩体结构面表面形态由宏观几何轮廓、表面形态和微观粗糙度三个要素组成，其中对结构面力学性质起决定性影响的是表面起伏形态。粗糙度是结构面起伏特征的定量描述，因此，通过结构面粗糙度来量化岩体结构表面。

目前方法存在的不足：

(1)粗糙度的参数统计和几何分形维数研究是建立在结构面二维剖面基础上提出的，在理论上存在一定程度的局限性，而且不同的研究者得到的结果差异较大，说明在二维面上进行粗糙度研究不能真实反映结构面的粗糙度系数。

(2)在三维结构面粗糙度的评价方面，研究理论较少，已有的研究结果还不能在工程中得到应用。

发明内容

针对现有技术存在的不足，本发明提出一种基于三角面积比法度量岩体结构面三维粗糙度的方法，以达到通过三维三角网格快速准确的测量结构面三维数据的目的。

一种基于三角面积比法度量岩体结构面三维粗糙度的方法，包括以下步骤：

步骤1、根据实际需求在被测岩体裸露面前放置一个范围杆采用摄相机在范围杆的一侧拍摄一张完整的岩体裸露面图片，在范围杆的另一侧拍摄另一张完整的岩体裸露面图片，并确定范围杆上圆盘与下圆盘之间的距离；

步骤2、采用罗盘测量被测岩体裸露面中延展范围大的结构面的倾向；

步骤3、采用三角面积比法获得被测岩体结构面的平面面积和表面积，具体如下：

步骤3-1、根据拍摄获得的两张岩体图片、范围杆上圆盘与下圆盘之间的距离和被测岩体裸露面中延展范围大的结构面的倾向，建立被测岩体的三维模拟模型；

步骤3-2、选定被测岩体的三维模拟模型的所有结构面，并确定所选结构面的边界点坐标和结构面出露部分与未出露部分交界点的坐标；

步骤3-3、采用最小二乘法线性回归方法对被测岩体结构面进行拟合，公式如下：

Z＝AX+BY+C(1)

其中，Z表示结构面竖向坐标；X表示结构面水平横向坐标；Y表示结构面水平纵向坐标；A、B和C表示方程系数；

步骤3-4、将所选的一个结构面的边界点坐标和交界点坐标代入步骤3-3拟合的公式中，获得公式如下：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i=nC+A\underset{i=1}{\overset{n}{\Σ}}{X}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{X}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i^2+B\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{Y}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right.---\left(2\right)$

其中，X_i为第i个点的水平横向坐标；Y_i为第i个点的水平纵向坐标；Z_i为第i个点的竖向坐标；n表示点的个数；

公式(2)的矩阵形式为：

$\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]\left[\begin{array}{c}C\\ A\\ B\end{array}\right]=\left[\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right]---\left(3\right)$

步骤3-5、判断公式(3)中 $\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]$ 矩阵的行列式是否为0，若是，则对应结构面与坐标轴垂直或平行，执行步骤3-7；否则，执行步骤3-6；

步骤3-6、根据克莱姆法则获得拟合公式的系数，并根据获得的系数计算求得对应结构面的倾向和倾角，公式如下：

$\mathrm{\α}={\tan }^{-1}\left(\sqrt{A^2+B^2}\right)---\left(4\right)$

其中，α表示倾向；

$\mathrm{\β}={\cos }^{-1}\left(\frac{A}{\sqrt{A^2+B^2}}\right)---\left(5\right)$

其中，β表示倾角；

步骤3-7、根据所选的一个结构面的边界点的X、Y坐标，获得该结构面在XY平面的投影面，并将该投影面分割成若干个三角形；

步骤3-8、确定每个三角形三个顶点的坐标，计算获得每个三角形的面积，并通过求和获得投影面的面积；

步骤3-9、根据该结构面的倾角和投影面的面积，获得结构面的平面面积；

步骤3-10、反复执行步骤3-4至步骤3-9，直至获得所有结构面的平面面积；

步骤4、将被测岩体的三维模拟模型中每个结构面进行三角形分割，并根据每个三角形的三个顶点坐标，计算每个三角形的面积，通过求和获得每个结构面的表面积；

计算每个三角形的面积公式如下：

$S_j=\frac{1}{2}\sqrt{I_1{I}_1+I_2{I}_2+I_3{I}_3}---\left(6\right)$

其中，S_j表示第j个三角形的面积；I₁＝(b₂-b₁)(c₃-c₁)-(c₂-c₁)(b₃-b₁)，I₂＝(c₂-c₁)(a₃-a₁)-(a₂-a₁)(c₃-c₁)，I₃＝(a₂-a₁)(b₃-b₁)-(b₂-b₁)(a₃-a₁)，其中，a₁、b₁、c₁分别表示三角形第一个顶点的坐标值，a₂、b₂、c₂分别表示三角形第二个顶点的坐标值，a₃，b₃，c₃分别表示三角形第三个顶点的坐标值；

步骤5、将每个结构面的平面面积除以结构面的表面积，获得结构面面积比；

步骤6、根据面积比获得结构面的粗糙度，公式如下：

JRC＝-4.6867SAR+6.9262(7)

其中，JRC表示结构面的粗糙度，SAR表示结构面面积比；

步骤7、根据被测岩体每个结构面的粗糙度，结合岩体被结构面切割的几何形态，确定该被测岩体中的关键岩块，以减少施工过程中工程事故出现的概率。

步骤3-6所述的根据克莱姆法则获得拟合公式的系数，公式如下：

C＝D₀/D，A＝D₁/D，B＝D₂/D(8)

其中，矩阵 $D=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_0=\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$

矩阵 $D_1=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_2=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right].$

本发明优点：

本发明一种基于三角面积比法度量岩体结构面三维粗糙度的方法，应用数字摄影测量系统进行结构面的采集，通过结构面出露部分的表面形态来反映整个结构面的三维空间形态，应用ShapeMetriX3D系统进行出露部分点数据的获取；通过最小二乘法对结构面进行平面拟合，计算结构面三角网格的总面积与平面面积之比，反应了拟合结构面的平整程度；该方法充分利用了摄影测量法可以快速准确的测量结构面三维数据的优势，在此基础上应用Delaunay三角剖分方法进行结构面表面的重建，实现了通过三维三角网格来反映结构面的表面形态；考虑二维粗糙度在表征结构面粗糙性的不足，在结构面表面三维重建的基础上，提出了比表面积法(surfacearearatio，SAR)，通过结构面的网格面积与拟合面积之间的比值来反映粗糙度大小。

附图说明

图1本发明一种实施例的基于三角面积比法度量岩体结构面三维粗糙度的方法流程图；

图2本发明一种实施例的二维平面分解图；

图3本发明一种实施例的二维平面分解三角形示意图；

图4本发明一种实施例的结构面平面面积的计算图；

图5本发明一种实施例的结构面投影图。

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。

一种基于三角面积比法度量岩体结构面三维粗糙度的方法，方法流程图如图1所示，包括以下步骤：

步骤1、根据实际需求在被测岩体裸露面前放置一个范围杆采用摄相机在范围杆的一侧拍摄一张完整的岩体裸露面图片，在范围杆的另一侧拍摄另一张完整的岩体裸露面图片，并确定范围杆上圆盘与下圆盘之间的距离；

本发明实施例中，所述的范围杆为ShapeMetrix3D三维模型重建软件系统附带产品，为同一厂家生产的配套装置，为现有公知产品。

步骤2、采用罗盘测量被测岩体裸露面中延展范围大的结构面的倾向；

步骤3、采用ShapeMetrix3D三维模型重建软件和三角面积比法，获得被测岩体结构面的平面面积和表面积，具体如下：

步骤3-1、将拍摄获得的两张岩体图片导入至三维模型重建软件中，并根据范围杆上圆盘与下圆盘之间的距离和被测岩体裸露面中延展范围大的结构面的倾向，建立被测岩体的三维模拟模型；

本发明实施例中，在现场岩体左右位置分别拍摄两张结构面照片，导入ShapeMetriX3D软件中，采用三维模型重建助手(SMXReconstructionAssistant)分别导入左右两张照片，初步合成岩体的3D图像，再采用SMXSurfaceTrimmer工具对边缘合成不佳的区域点进行编辑删除，获得直观、清晰、有效的出露结构面，即获得岩体表面的三维立体图像(被测岩体的三维模拟模型)。

获得的三维立体图像与实际岩体的大小不同，因此需要将三维图像真实化，具体为：打开图像，在图片上找到圆盘的中心点，根据步骤1测量的这两点间标杆的真实距离，并将测量的两圆盘间距离输入至三维立体图像中，则此时三维图上任意两点的距离已经与实际相吻合；并采用罗盘量测一个区域大的结构面，在三维立体图像中圈出对应的结构面，并输入该面的倾向，将结构面的倾角真实化。

步骤3-2、选定被测岩体的三维模拟模型的所有结构面，并确定所选结构面的边界点坐标和结构面出露部分与未出露部分交界点的坐标；

本发明实施例中，对获得的岩体三维立体图像进行结构面信息的分析(JMXAnalyst)，具体为：在图像上选择同一个结构面上的点进行标记，ShapeMetriX3D系统自动实现每个结构面的识别、定位以及几何形态信息真实参数(产状、迹长、间距、断距)的获取。

为了获得结构面的点群数据，选择延展大的结构面，手动选择结构面出露部分的点，导出这些点的三维坐标，选择点遵循两条规则：1、选择的点包括边界上的点以及出露与未出露部分交界上的点；2、选择点均匀，包括起伏部分等特殊位置的点。选择完成后，采用ShapeMetriX3D系统以excel形式保存选择的点群的三维坐标。

步骤3-3、采用最小二乘法线性回归方法对被测岩体结构面进行拟合；

平面形态拟合方法有最小二乘法和特征值法，本发明实施例中，应用最小二乘法线性回归的方法进行结构面拟合。在三维直角坐标系下，设结构面的平面方程为：

Z＝AX+BY+C(1)

其中，Z表示结构面竖向坐标；X表示结构面水平横向坐标；Y表示结构面水平纵向坐标；A、B和C表示方程系数；

步骤3-4、将所选的一个结构面的边界点坐标和交界点坐标代入步骤3-3拟合的公式中；

本发明实施例中，通过ShapeMetriX3D系统获得单一结构面的点群数据后，对结构面进行平面拟合，计算出结构面的倾向，并与ShapeMetriX3D所得到的结构面倾向进行比较，对结构面的点选择进行校核，然后计算结构面平面的面积。

通过摄影测量方法获得结构面上不共线的n(n≥3)个点的物方空间坐标为：

P_i(X_i，Y_i，Z_i)，i＝1，2，……n，

其中，P_i表示第i个点；X_i为第i个点的水平横向坐标；Y_i为第i个点的水平纵向坐标；Z_i为第i个点的竖向坐标；n表示点的个数；

采用线性回归的方法求解方程系数A，B，C，公式为：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i=nC+A\underset{i=1}{\overset{n}{\Σ}}{X}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{X}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i^2+B\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{Y}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right.---\left(2\right)$

公式(2)的矩阵形式为：

$\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]\left[\begin{array}{c}C\\ A\\ B\end{array}\right]=\left[\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right]---\left(3\right)$

步骤3-5、判断公式(3)中 $\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]$ 矩阵的行列式是否为0，若是，则对应结构面与坐标轴垂直(此时，倾角为90°)或平行(此时，倾角为0°)，执行步骤3-7；否则，执行步骤3-6；

步骤3-6、根据克莱姆法则获得拟合公式的系数，并根据获得的系数计算求得对应结构面的倾向和倾角；

根据克莱姆法则获得拟合公式的系数，公式如下：

C＝D₀/D，A＝D₁/D，B＝D₂/D(8)

其中，矩阵 $D=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_0=\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$

矩阵 $D_1=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_2=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right]$

根据平面的方程，结构面的法线向量为(A，B，-1)，设结构面的倾向为β，倾角为α，则根据结构面的法向量按照下式计算倾向和倾角：

$a={\tan }^{-1}\left(\sqrt{A^2+B^2}\right)---\left(4\right)$

其中，α表示倾向；

$\mathrm{\β}={\cos }^{-1}\left(\frac{A}{\sqrt{A^2+B^2}}\right)---\left(5\right)$

其中，β表示倾角；

步骤3-7、根据所选的一个结构面的边界点的X、Y坐标，获得该结构面在XY平面的投影面，并将该投影面分割成若干个三角形；

本发明实施例中，对投影面进行Delaunay三角剖分，所述的Delaunay三角剖分方法是将点群坐标投影到二维平面上，应用三角生长法进行结构面构建，获得三角形网格的点、边、面之间的拓扑关系，并转化到三维点群上，实现基于摄影测量结构面信息，结构面表面形态的三维网格模拟。

步骤3-8、确定每个三角形三个顶点的坐标，计算获得每个三角形的面积，并通过求和获得投影面的面积；

本发明实施例中，如图2所示，在二维平面上结构面的投影为多边形，通过计算每个三角形的面积来求出总的面积；如图3所示，每个三角形的三个顶点坐标为M(x₁，y₁)，N(x₂，y₂)，Q(x₃，y₃)，则三角形的面积S_j′为：

$S_{j^{\′}}=\frac{1}{2}\left[\begin{array}{ccc}{x}_1& y_1& 1\\ x_2& y_2& 1\\ x_3& y_3& 1\end{array}\right]---\left(9\right)$

则结构面在二维平面上的面积(投影面积)S_p1的面积计算公式为：

$S_{p1}=\underset{j^{\′}=1}{\overset{m}{\Σ}}{S}_{j^{\′}}---\left(10\right)$

其中，S_p1表示结构面在二维平面(投影面)上的面积，S_j′表示二维平面(投影面)上三角网格的面积，m表示结构面在二维平面三角网格的个数。

步骤3-9、根据该结构面的倾角和投影面的面积，获得结构面的平面面积；

本发明实施例中，如图4所示，根据该结构面的倾角α和投影面的面积S_p1，获得结构面的平面面积S_p，公式为：

步骤3-10、反复执行步骤3-4至步骤3-9，直至获得所有结构面的平面面积；

步骤4、将被测岩体的三维模拟模型中每个结构面进行Delaunay三角剖分，并根据每个三角形的三个顶点坐标，计算每个三角形的面积，通过求和获得每个结构面的表面积；

本发明实施例中，结构面通过三角剖分之后，表面形态由k个三角形网格构成，每个三角形网格对应的点坐标具有拓扑关系，因此，结构面表面的面积(S_t)通过计算每个三维三角网格的面积然后求和即可得到。其中，一个三角网格的三个顶点坐标分别为：(a₁，b₁，c₁)，(a₂，b₂，c₂)，(a₃，b₃，c₃)，则三角网格的面积S_j为：

$S_j=\frac{1}{2}\sqrt{I_1{I}_1+I_2{I}_2+I_3{I}_3}---\left(6\right)$

其中，S_j表示第j个三角形的面积；I₁＝(b₂-b₁)(c₃-c₁)-(c_2-c₁)(b₃-b₁)，I₂＝(c₂-c₁)(a₃-a₁)-(a₂-a₁)(c₃-c₁)，I₃＝(a₂-a₁)(b₃-b₁)-(b₂-b₁)(a₃-a₁)，其中，a₁、b₁、c₁分别表示三角形第一个顶点的坐标值，a₂、b₂、c₂分别表示三角形第二个顶点的坐标值，a₃，b₃，c₃分别表示三角形第三个顶点的坐标值；

结构面表面的面积S_t的计算公式为：

步骤5、将每个结构面的平面面积除以结构面的表面积，获得结构面面积比；

本发明实施例中，如图5所示，在XYZ坐标系中，上端具有分割三角形的图形为结构面，下端为结构面的平面面积；三角面积比法(surfacearearatio，SAR)，结合结构面的起伏形态，应用结构面的平面面积与起伏三角形的面积和的比来表征三维结构面的粗糙度，表面积比公式如下：

$SAR=\frac{S_p}{S_t}\×100\%---\left(11\right)$

步骤6、根据面积比获得结构面的粗糙度；

本发明实施例中，根据现场结构面的测量，应用修正直边法进行了粗糙度系数的计算，对SAR和粗糙度系数JRC之间的关系，做了相关性分析，根据现场的测量得出了JRC与SAR的关系，即：

JRC＝-4.6867SAR+6.9262(7)

其中，JRC表示结构面的粗糙度，SAR表示结构面面积比；

步骤7、根据被测岩体每个结构面的粗糙度，结合岩体被结构面切割的几何形态，确定该被测岩体中的关键岩块，以减少施工过程中工程事故出现的概率。

本发明实施例中，采用GeoSMA-3D关键块体三维分析系统，在GeoSMA-3D系统中采用圆盘模型进行结构面的模拟，将岩体被结构面切割的几何形态参数(中心点坐标、倾向、倾角、圆盘半径、内摩擦角和黏聚力)按照顺序存储在txt文档中，程序从txt文档中调用数据进行结构面的模拟；利用上述结构面粗糙度计算值，结合岩土工程建模与分析系统软件，对实际工程的关键块体进行识别，找出工程中的关键块，从而在施工过程中对关键块部位进行特殊处理，尽量减少工程事故出现的概率。

Claims (2)

1.一种基于三角面积比法度量岩体结构面三维粗糙度的方法，其特征在于，包括以下步骤：

步骤1、根据实际需求在被测岩体裸露面前放置一个范围杆，采用摄相机在范围杆的一侧拍摄一张完整的岩体裸露面图片，在范围杆的另一侧拍摄另一张完整的岩体裸露面图片，并确定范围杆上圆盘与下圆盘之间的距离；

步骤2、采用罗盘测量被测岩体裸露面中延展范围大的结构面的倾向；

步骤3、采用三角面积比法获得被测岩体结构面的平面面积和表面积，具体如下：

步骤3-1、根据拍摄获得的两张岩体图片、范围杆上圆盘与下圆盘之间的距离和被测岩体裸露面中延展范围大的结构面的倾向，建立被测岩体的三维模拟模型；

步骤3-2、选定被测岩体的三维模拟模型的所有结构面，并确定所选结构面的边界点的坐标和结构面出露部分与未出露部分交界点的坐标；

步骤3-3、采用最小二乘法线性回归方法对被测岩体结构面进行拟合，公式如下：

Z＝AX+BY+C(1)

其中，Z表示结构面竖向坐标；X表示结构面水平横向坐标；Y表示结构面水平纵向坐标；A、B和C表示方程系数；

步骤3-4、将所选的一个结构面的边界点坐标和交界点坐标代入步骤3-3拟合的公式中，获得公式如下：

$\left\{\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i=nC+A\underset{i=1}{\overset{n}{\Σ}}{X}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{X}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i^2+B\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i=C\underset{i=1}{\overset{n}{\Σ}}{Y}_i+A\underset{i=1}{\overset{n}{\Σ}}{X}_i{Y}_i+B\underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right.---\left(2\right)$

其中，X_i为第i个点的水平横向坐标；Y_i为第i个点的水平纵向坐标；Z_i为第i个点的竖向坐标；n表示点的个数；

公式(2)的矩阵形式为：

$\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]\left[\begin{array}{c}C\\ A\\ B\end{array}\right]=\left[\begin{array}{c}\underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right]---\left(3\right)$

步骤3-5、判断公式(3)中 $\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right]$ 矩阵的行列式是否为0，若是，则对应结构面与坐标轴垂直或平行，执行步骤3-7；否则，执行步骤3-6；

步骤3-6、根据克莱姆法则获得拟合公式的系数，并根据获得的系数计算求得对应结构面的倾向和倾角，公式如下：

$\mathrm{\α}={\tan }^{-1}\left(\sqrt{A^2+B^2}\right)---\left(4\right)$

其中，α表示倾向；

$\mathrm{\β}={\cos }^{-1}\left(\frac{A}{\sqrt{A^2+B^2}}\right)---\left(5\right)$

其中，β表示倾角；

步骤3-7、根据所选的一个结构面的边界点的X、Y坐标，获得该结构面在XY平面的投影面，并将该投影面分割成若干个三角形；

步骤3-8、确定每个三角形三个顶点的坐标，计算获得每个三角形的面积，并通过求和获得投影面的面积；

步骤3-9、根据该结构面的倾角和投影面的面积，获得结构面的平面面积；

步骤3-10、反复执行步骤3-4至步骤3-9，直至获得所有结构面的平面面积；

步骤4、将被测岩体的三维模拟模型中每个结构面进行三角形分割，并根据每个三角形的三个顶点坐标，计算每个三角形的面积，通过求和获得每个结构面的表面积；

计算每个三角形的面积公式如下：

$S_j=\frac{1}{2}\sqrt{I_1{I}_1+I_2{I}_2+I_3{I}_3}---\left(6\right)$

其中，S_j表示第j个三角形的面积；I₁＝(b₂-b₁)(c₃-c₁)-(c₂-c₁)(b₃-b₁)，I₂＝(c₂-c₁)(a₃-a₁)-(a₂-a₁)(c₃-c₁)，I₃＝(a₂-a₁)(b₃-b₁)-(b₂-b₁)(a₃-a₁)，其中，a₁、b₁、c₁分别表示三角形第一个顶点的坐标值，a₂、b₂、c₂分别表示三角形第二个顶点的坐标值，a₃，b₃，c₃分别表示三角形第三个顶点的坐标值；

步骤5、将每个结构面的平面面积除以结构面的表面积，获得结构面面积比；

步骤6、根据面积比获得结构面的粗糙度，公式如下：

JRC＝-4.6867SAR+6.9262(7)

其中，JRC表示结构面的粗糙度，SAR表示结构面面积比；

步骤7、根据被测岩体每个结构面的粗糙度，结合岩体被结构面切割的几何形态，确定该被测岩体中的关键岩块。

2.根据权利要求1所述的基于三角面积比法度量岩体结构面三维粗糙度的方法，其特征在于，步骤3-6所述的根据克莱姆法则获得拟合公式的系数，公式如下：

C＝D₀/D，A＝D₁/D，B＝D₂/D(8)

其中，矩阵 $D=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_0=\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$

矩阵 $D_1=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i^2\end{array}\right],$ 矩阵 $D_2=\left[\begin{array}{ccc}n& \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{X}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i^2& \underset{i=1}{\overset{n}{\Σ}}{X}_i{Z}_i\\ \underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{X}_i\underset{i=1}{\overset{n}{\Σ}}{Y}_i& \underset{i=1}{\overset{n}{\Σ}}{Y}_i{Z}_i\end{array}\right].$