一种由单节理迹线生成曲面节理圆盘的岩体结构重构方法

1.本发明涉及岩体离散裂隙网络领域，尤其涉及一种由单节理迹线生成曲面节理圆盘的岩体结构重构方法。


背景技术：

2.岩体由节理和完整的岩块组成，节理是指在岩体内部发育的具有一定方向、规模和形态的物质分界面或不连续面，由于节理的存在，使岩体的力学和水力特性非常复杂。在自然界中，岩石材料具有遮光特性，岩体内部节理面的几何特征难以测得，只能通过岩体露头与节理面的交线，即迹线来推测。因此为了研究岩体的稳定以及渗流等特性，通过迹线信息构造岩体的三维离散裂隙网络模型(dfn)是最有效的方法。
3.圆盘节理模型因其简单的参数和良好的模拟结果而成为应用最广泛的dfn之一，此模型将岩体节理面假设为圆盘形状，由节理参数模拟生成一个节理圆盘，重复单个节理圆盘的生成过程而构成节理圆盘模型。然而节理圆盘的生成大多假设节理是完全平面，这与实际情况不符。自然形成的节理面往往是具有粗糙起伏的不规则几何面，其表面形态非常复杂，具有很强的随机性，综上所述，现有方法中构建的平面节理圆盘无法准确模拟真实的岩体结构，导致测量的岩体特性不准确。


技术实现要素：

4.本发明提供一种由单节理迹线生成曲面节理圆盘的岩体结构重构方法，以克服现有方法中构建的平面节理圆盘无法准确模拟真实的岩体结构，测量岩体特性不准确的问题。
5.为了实现上述目的，本发明的技术方案是：
6.一种由单节理迹线生成曲面节理圆盘的岩体结构重构方法，如图1所示，包括以下步骤：
7.s1、基于三维点云技术提取岩体的节理信息，计算节理的产状；
8.s2、根据所述节理信息确定节理的二维剖面线，将所述节理二维剖面线的中心线作为x轴，将所述中心线的垂线作为z轴，以所述节理二维剖面线的端点在所述x轴上的投影为原点，在垂直面上建立局部直角坐标系，并计算所述节理二维剖面线的粗糙度参数，所述粗糙度参数包括rq和hurst指数；
9.s3、根据所述粗糙度参数，采用矩阵变换的方法，生成一个曲面节理圆盘；
10.s4、将所述曲面节理圆盘进行坐标转换，包括移动到节理面位置和旋转至与迹线拟合节理面产状相同的位置，使所述曲面节理圆盘的坐标与全局坐标系一致；
11.s5、根据坐标转换后的曲面节理圆盘，使用三维节理面网络模拟方法构建岩体的非平面节理圆盘模型，以研究岩体的稳定以及渗流特性，并对岩体结构进行评价。
12.进一步地，所述步骤s3包括：
13.s31、根据hurst指数计算特征值l，并分别计算rho矩阵和phi矩阵，根据节理面的
精度常数n生成行矩阵x和行矩阵y；
14.s32、生成一个与所述rho矩阵和phi矩阵大小相同的随机矩阵c，对所述随机矩阵c进行二维离散傅里叶变换得到矩阵g；
15.s33、将矩阵g做重新排列傅里叶变换得到矩阵g
′
，所述矩阵g
′
与所述phi矩阵做比值运算得到矩阵w，将所述矩阵w做重新排列傅里叶变换得到矩阵w
′
；
16.s34、对所述矩阵w
′
做二维离散傅里叶逆变换得到矩阵k；
17.s35、根据所述矩阵k中每个元素的实部得到矩阵t，根据所述矩阵t计算出分形维数矩阵；
18.s36、计算所述分形维数矩阵内元素的标准差μ，并根据所述标准差μ得到矩阵z，根据所述行矩阵x、行矩阵y和矩阵z中的坐标点数据建立方形曲面，所述方形曲面的内切圆为曲面节理圆盘。
19.进一步地，所述步骤s1包括以下步骤：
20.s11、基于数字近景摄影测量技术建立岩体的纹理化三维点云模型；
21.s12、提取一条节理迹线作为节理三维轮廓线，采用最小二乘法将所述节理迹线的拐点坐标拟合得到拟合平面，将所述拟合平面作为节理面，并计算所述节理面的法向量；
22.s13、根据所述节理面的法向量计算该节理的产状。
23.进一步地，所述步骤s13中的产状包括倾向和倾角；
24.所述节理迹线的拐点所拟合的节理面的法向量为n，所述倾向为y轴与所述法向量n在水平投影之间的夹角，范围为0-360
°
；所述倾角为水平面与所述节理面之间的二面角，范围为0-90
°
；所述y轴为地球坐标系下的正北方向。
25.进一步地，所述步骤s2并计算所述节理二维剖面线的rq和hurst指数中粗糙度参数rq是所述节理二维剖面线上点的高度z相对于x轴的均方根，rq的计算公式为：
[0026][0027]
公式(1)中：z为节理二维剖面线与x轴之间的距离，l为节理二维剖面线在x轴方向上的最大投影长度；
[0028]
所述步骤s2中hurst指数采用盒子分形维数的方法计算，包括以下步骤：
[0029]
s21、将节理二维剖面线置于采样窗口内，使用固定比率成几何级数的不同尺寸盒子将采样窗口离散成多个方形盒子，盒子的尺寸r分别为：l/128、l/64、l/32、l/16和l/8；
[0030]
s22、通过计算完全覆盖采样窗口的节理二维剖面线的盒子的数量，得出盒子数量s与盒子尺寸倒数1/r在双对数坐标图中的关系曲线；
[0031]
s23、采用最小二乘法对所述关系曲线进行拟合，得到拟合直线，确定拟合直线的斜率d；
[0032]
s24、根据h＝2-d确定hurst指数，其中h为hurst指数，d为拟合直线的斜率。
[0033]
进一步地，所述步骤s4包括以下步骤：
[0034]
s41、对所述曲面节理圆盘上的点进行平移，将圆心平移至迹线的中点，平移公式为：(x，y，z)＝(x+x0，y+y0，z+z0)，其中(x，y，z)为生成曲面节理圆盘的点坐标，(x0，y0，z0)为全局坐标系下该迹线上所有拐点的平均坐标值；
[0035]
s42、建立如下坐标转换矩阵：
[0036][0037]
公式(2)中：γ＝dipx-α，dipx是全局坐标系x轴方向上的倾向，α是节理面的倾向，β是节理面的倾角；
[0038]
所述坐标转换矩阵nto的逆矩阵otn为：
[0039][0040]
将曲面节理圆盘上的点坐标乘以坐标转换矩阵的逆矩阵otn，将节理圆盘旋转至与迹线拟合节理面产状相同的位置。有益效果：本发明基于近景摄影测量建立的现场岩体三维点云模型，通过提取迹线数据获取岩体的节理信息，用分形方法重建节理表面的粗糙度，提出了一种由单条节理迹线生成一个曲面节理圆盘的方法，相比于传统节理圆盘，构建曲面节理圆盘考虑了节理粗糙度的因素，增加了节理圆盘的真实性和适用性，可以更加准确的模拟真实的岩体结构，测得的岩体特性更准确。
附图说明
[0041]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
[0042]
图1为本发明生成曲面节理圆盘的方法流程图；
[0043]
图2为三维点云模型中提取一条迹线的示意图；
[0044]
图3为选取一条迹线准备处理的示意图；
[0045]
图4为通过提取迹线获取其节理几何参数的示意图；
[0046]
图5为计算节理产状的示意图；
[0047]
图6为节理二维剖面线示意图；
[0048]
图7为盒子大小取所有值时的盒子网络示意图；
[0049]
图8为盒子数量与盒子尺寸倒数在双对数坐标图中的关系图；
[0050]
图9为生成一个曲面节理圆盘的示意图；
[0051]
图10为坐标转换后的曲面节理圆盘示意图；
[0052]
图11为非平面节理圆盘模型示意图。
具体实施方式
[0053]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0054]
本实施例提供了一种由单节理迹线生成曲面节理圆盘的岩体结构重构方法，包括
以下步骤：
[0055]
s1、基于三维点云技术提取岩体的节理信息，计算节理的产状；
[0056]
具体地，步骤s1包括如下步骤：
[0057]
s11、基于数字近景摄影测量技术建立岩体的纹理化三维点云模型；
[0058]
具体地，通过数字近景摄影测量技术建立实际边坡岩体的纹理化岩体三维点云模型，在岩体露头上，节理面与岩石的表面交切形成迹线，节理的信息通过提取迹线的信息获得。由于节理是具有粗糙起伏的不规则几何曲面，所以迹线是一条不规则的曲线。在本实施例中，在地球坐标系下，以正东为x轴，正北为y轴，重力的反方向为z轴，建立全局空间坐标系，提取并记录模型内迹线上拐点的坐标，为充分描述迹线的空间特征，应提取尽可能多的拐点，如图2所示为三维点云模型中提取出一条迹线的示意图，图2中a为提取出的一条迹线。
[0059]
s12、提取一条节理迹线作为节理三维轮廓线，采用最小二乘法将所述节理迹线的拐点坐标拟合得到拟合平面，将所述拟合平面作为节理面；
[0060]
具体地，选取一条节理迹线，将迹线作为节理三维轮廓线，如图3所示为选取的一条迹线示意图，如图4所示为一条迹线获取节理几何参数的示意图，图4中b为从3d模型中提取出的一条迹线，c为与拟合平面垂直的投影面，d为拟合平面，e为垂直面上的投影线，根据迹线拐点的坐标采用最小二乘法拟合得到拟合平面，拟合平面作为迹线的节理面，根据节理面的法向量计算该节理的产状，即倾向和倾角。设定一个与节理面垂直的平面，将节理三维轮廓线在垂直平面中形成投影线，移动平面使垂直平面中投影线的长度达到最大，将此时的投影线作为节理二维剖面线。
[0061]
s13、根据所述节理面的法向量计算该节理的产状；
[0062]
具体地，如图5所示为节理产状计算示意图，产状的计算方法为，拐点p1、p2和p3所拟合的节理面对应图5中虚线圆圈所示平面，其节理面的法向量为n(xn，yn，zn)，拟合平面的产状即倾向和倾角，其中倾向定义为正北方向与法向量n在水平投影之间的夹角，范围为0-360
°
，对应图5中的β；倾角定义为水平面与节理平面之间的二面角，范围为0-90
°
，对应图5中的α，倾向和倾角可根据法向量n进行推断。
[0063]
s2、根据所述节理信息确定节理的二维剖面线，将所述节理二维剖面线的中心线作为x轴，将所述中心线的垂线作为z轴，以所述节理二维剖面线的端点在所述x轴上的投影为原点，在垂直面上建立局部直角坐标系，并计算所述节理二维剖面线的粗糙度参数，所述粗糙度参数包括rq和hurst指数；
[0064]
具体地，步骤s2中rq是节理二维剖面线上点的高度z0相对于x轴的均方根，其计算公式为：
[0065][0066]
公式(1)中：z0为节理二维剖面线上点的高度，即该点与x轴之间的距离，l为节理二维剖面线在x轴方向上的最大投影长度，如图6所示为节理二维剖面线示意图，图6中节理二维剖面线的粗糙度参数rq＝0.04m，l＝7.75m。
[0067]
步骤s2中hurst指数采用盒子分形维数的方法进行计算，具体包括以下步骤：
[0068]
s21、将节理二维剖面线置于采样窗口内，使用固定比率成几何级数的不同尺寸盒
子将采样窗口离散成多个方形盒子，盒子的尺寸r分别为：l/128、l/64、l/32、l/16、l/8；
[0069]
s22、通过计算完全覆盖采样窗口的节理二维剖面线的盒子的数量，得出盒子数量s与盒子尺寸倒数1/r在双对数坐标图中的关系曲线；
[0070]
s23、采用最小二乘法对所述关系曲线进行拟合，得到拟合直线，确定拟合直线的斜率d；
[0071]
s24、根据公式：h＝2-d确定hurst指数，其中h为hurst指数，d为拟合直线的斜率。
[0072]
具体地，采用盒子分形维数的方法计算hurst指数：首先将节理二维剖面线置于采样窗口内，使用不同尺寸大小的盒子将采样窗口离散成小的方形盒子。其中采样窗口中盒子数量为整数，盒子尺寸为固定比率的几何级数，所取得盒子尺寸r分别为l/128，l/64，l/32，l/16，l/8，其中l为节理二维剖面线在x轴方向最大投影长度。如图7所示为盒子大小取所有值时的盒子网络，当盒子大小取不同值时，计算完全覆盖采样窗口中节理二维剖面线的盒子的数量，并绘制盒子数量s与盒子尺寸倒数1/r在双对数坐标图中的关系曲线，将曲线用最小二乘法计算得到拟合直线，拟合直线的斜率记为d。如图8所示为s和1/r在双对数坐标图中的关系。图中所示的d＝1.156，根据d与h的关系：h＝2-d，从而求得h，h＝0.084。
[0073]
s3、生成一个曲面节理圆盘；
[0074]
具体地，步骤s3包括以下步骤：
[0075]
s31、根据hurst指数计算特征值l，并分别计算rho矩阵和phi矩阵，根据节理面的精度常数n生成行矩阵x和行矩阵y；
[0076]
具体地，做初始数据处理，首先计算特征值l，l为hurst指数的倒数，即l＝1/h，然后生成行矩阵x和y：
[0077]
x＝y＝[-1 2-n-1 2
1-n-1 2
2-n-1 ... 1]；
ꢀꢀꢀ
(2)
[0078]
矩阵x和矩阵y的大小为1
×
(2
n+1
+1)，是元素从-1到1，公差为2-n
的等差数列，n为节理面精度常数，n越大精度越大，通常为8-11；将矩阵x分块得到分块矩阵xx：
[0079][0080]
分块矩阵xx的子矩阵为矩阵x，大小为(2
n+1
+1)
×
1，然后计算rho矩阵和phi矩阵：
[0081][0082]
phi＝rho
(1+l)
；
ꢀꢀꢀ
(5)
[0083]
公式(4)和(5)中：abs函数表示将矩阵内的元素进行绝对值计算的函数，rho矩阵和phi矩阵的大小为(2
n+1
+1)
×
(2
n+1
+1)。
[0084]
s32、生成一个与所述rho矩阵和phi矩阵大小相同的随机矩阵c，对所述随机矩阵c进行二维离散傅里叶变换得到矩阵g；
[0085]
具体地，生成一个大小为(2
n+1
+1)
×
(2
n+1
+1)的随机矩阵c，其中matlab内置函数为fft2，对随机矩阵c进行二维离散傅里叶变换得到矩阵g，矩阵g内各个元素的计算公式为：
[0086][0087]
[0088]
公式(6)和(7)中：g
p+1，q+1
为矩阵g的第p+1行第q+1列元素，c
j+1，k+1
为矩阵c的第j+1行第k+1列元素，是单位复根，i是虚数单位，p、q、j、k的范围为0到2
n+1
。
[0089]
s33、将矩阵g做重新排列傅里叶变换得到矩阵g
′
，所述矩阵g
′
与所述phi矩阵做比值运算得到矩阵w，将所述矩阵w做重新排列傅里叶变换得到矩阵w
′
；
[0090]
具体地，首先将矩阵g写为分块矩阵形式：
[0091][0092]
公式(8)中：子矩阵g1、g2、g3、g4的大小分别为2n×2n
、(2n)
×
(2n+1)、(2n+1)
×
(2n)、(2n+1)
×
(2n+1)；
[0093]
将矩阵g进行重新排列傅里叶变换得到矩阵g
′
，其中matlab内置函数为fftshift，矩阵g
′
表示为：
[0094][0095]
继续得到矩阵w，并将矩阵w写成分块矩阵的形式，表示为：
[0096][0097]
公式(10)中：矩阵w的大小为(2
n+1
+1)
×
(2
n+1
+1)，四个子矩阵w1、w2、w3、w4的大小分别为2n×2n
、(2n)
×
(2n+1)、(2n+1)
×
(2n)、(2n+1)
×
(2n+1)；
[0098]
接着将矩阵w进行重新逆排列傅里叶变换得到矩阵w
′
，其中matlab内置函数为iffshift，矩阵w
′
表示为：
[0099][0100]
s34、对所述矩阵w
′
做二维离散傅里叶逆变换得到矩阵k；
[0101]
具体地，对大小为(2
n+1
+1)
×
(2
n+1
+1)的矩阵w
′
进行二维离散傅里叶逆变换，其中matlab内置函数为ifft2，得到大小为(2
n+1
+1)
×
(2
n+1
+1)的矩阵k，矩阵k内各个元素的计算公式为：
[0102][0103][0104]
公式(12)和(13)中：k
p，q
为矩阵k的第p行第q列元素，w
′
j，k
为矩阵w的第j行第k列元素，是单位复根，i是虚数单位，p、q、j、k的范围为1到2
n+1
+1。
[0105]
s35、根据所述矩阵k中每个元素的实部得到矩阵t，根据所述矩阵t计算出分形维数矩阵roughsurf；
[0106]
具体地，只记录矩阵k中每个元素的实部得到大小为(2
n+1
+1)
×
(2
n+1
+1)的矩阵t，即t＝real(k)，其中real函数表示返回矩阵k中每个元素的实部；由矩阵t计算得到大小为(2
n+1
+1)
×
(2
n+1
+1)的分形维数矩阵roughsurf：
[0107]
roughsurf＝t-t
2n+1，2n+1
。
ꢀꢀꢀ
(14)
[0108]
s36、计算所述分形维数矩阵roughsurf内元素的标准差μ，并根据所述标准差μ得到矩阵z，根据所述行矩阵x、行矩阵y和矩阵z中的坐标点数据建立方形曲面，以所述方形曲面的坐标原点为圆心、所述节理二维剖面线在x轴上的最大投影距离为直径，建立曲面节理圆盘；
[0109]
具体地，分形维数矩阵roughsurf内元素的标准差μ的计算公式为：
[0110][0111]
公式(15)中：i、j的范围为1到2
n+1
+1；
[0112]
大小为(2
n+1
+1)
×
(2
n+1
+1)的矩阵z的计算公式为：
[0113][0114]
公式(16)中：rq是节理二维剖面线上点的高度z相对于x轴的均方根；
[0115]
具体地，将行矩阵x、行矩阵y内的元素扩大l/2倍，得到矩阵x
′
和y
′
，其中l表示所述节理二维剖面线在x轴方向上的最大投影长度，即：
[0116][0117]
具体地，将矩阵z内(2
n+1
+1)
×
(2
n+1
+1)个元素作为相应数量的点的z坐标值，将矩阵x
′
第n列元素值作为矩阵z第n列元素对应点的x坐标，将矩阵y
′
第n列的元素值作为矩阵z第n行元素对应点的y坐标，即可得到(x，y，z)的坐标点数据，接着由坐标点数据建立方形曲面，如图9所示，该方形曲面以坐标原点为圆心，以l为直径的圆即为曲面节理圆盘。
[0118]
s4、将所述曲面节理圆盘进行坐标转换，包括移动到节理面位置和旋转至与迹线拟合节理面产状相同的位置，使所述曲面节理圆盘的坐标与全局坐标系一致；
[0119]
具体地，步骤s4包括以下步骤：
[0120]
s41、平移；对所述曲面节理圆盘上的点进行平移，将圆心平移至迹线的中点，平移公式为：(x，y，z)＝(x+x0，y+y0，z+z0)，其中(x，y，z)为生成曲面节理圆盘的点坐标，(x0，y0，z0)为全局坐标系下该迹线上所有拐点的平均坐标值；
[0121]
s42、旋转；建立如下坐标转换矩阵：
[0122][0123]
其中：γ＝dipx-α，dipx是全局坐标系x轴方向上的倾向，α是节理面的倾向，β是节理面的倾角；
[0124]
所述坐标转换矩阵nto的逆矩阵otn为：
[0125][0126]
将曲面节理圆盘上的点坐标乘以坐标转换矩阵的逆矩阵otn，将节理圆盘旋转至与迹线拟合节理面产状相同的位置，如图10所示为坐标转换后的曲面节理圆盘。
[0127]
s5、根据坐标转换后的曲面节理圆盘，使用三维节理面网络模拟的方法构建岩体
的非平面节理圆盘模型，如图11所示，将岩体结构进行仿真重构后即可研究岩体的稳定性以及渗流等特性，并可以依此仿真数据对岩体结构进行进一步评价，以便后续进行其他研究。
[0128]
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。