一种基于分形算法的区域地质体岩体裂隙模拟方法与流程

本发明涉及岩体裂隙模拟方法，尤其是隧道工程区裂隙模拟方法，具体为一种基于分形算法的区域地质体岩体裂隙模拟方法。

背景技术：

岩体中裂隙是地质体中常见的地质构造，其形态复杂而且千奇百怪。用计算机模拟它们生成的时候，其建模是很困难的。为了能够生成形态逼真的裂隙模型，一般来说，用简单的递归方法是不能从根本上解决这个问题的。分形方法的引入，为裂隙模型构型提供了新的途径。随机LS文法引入了随机变量使表象不同形态的规则按不同的概率分布随机出现，从而产生更加自然、生动的图形来。然而在传统的LS文法中，只能掌握其大概的形成规则，在规定了生长方式之后，只能按照随机的方法来绘制分形裂隙，而对于一些具体的细节却无法自由调节。譬如:各种生长方方式的比例，以及分形元的比例。目前的分形研究局限于二维平面，实际上裂隙分布在三维空间，因此现在亟待研究一种方法真实模拟三维岩体裂隙。

技术实现要素：

本发明所要解决的技术问题是，克服现有技术的缺点，提供一种基于分形算法的区域地质体岩体裂隙模拟方法，本发明设计的岩体裂隙模拟方法是根据随机LS文法生成，提出了参数化设计的思想，研究了如何在生成分形裂隙时，利用生长角度参数a和概率控制参数U、V来控制裂隙树形，可以绘制形态各异的三维裂隙分形裂隙；本发明指出利用这些参数可以更好地调节各种分形元的比例，从而达到人为控制三维裂隙分形裂隙树形的目的，并通过程序实现了LS文法绘制三维裂隙分形裂隙的参数化设计，仿真模拟了该隧道工程区的三维裂隙。

本发明的目的是通过隧道工程区的参数计算，模拟真实空间分布工程区三维裂隙信息。

本发明的应用基础是区域地质体三维模型，该地质体三维模型分区块包括三维坐标存储、隧道围岩属性信息，以及断层、产状等区域地质调查数据。

本发明提供了一种基于分形算法的区域地质体岩体裂隙模拟方法，包含以下步骤：

步骤一、根据隧道围岩属性划分不同区域：

在隧道路线设计图中，含隧道起始点里程L₀，坐标(x₀，y₀,z₀)，隧道终止里程L_n，坐标(x_n，y_n,z_n)，中间部分隧道路线各节点里程值为L₁，L₂，L₃，…，L_n-1，坐标为(x₁，y₁,z₁)，(x₂，y₂,z₂)，(x₃，y₃,z₃)，…，(x_n-1，y_n-1,z_n-1)。隧道整体区域以隧道不同的围岩属性信息分为不同的里程信息l₀，l₂，l₃，…，l_m首先根据里程信息计算不同里程分界点的三维坐标(x_i，y_i,z_i)(0<i<m)：

x_i＝x_j+(x_j+1-x_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)，

y_i＝y_j+(y_j+1-y_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)，

z_i＝z_j+(z_j+1-z_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)；

步骤二、根据不同区域的围岩属性信息计算此区域分形深度：

隧道围岩属性根据破碎程度及岩性分为5级，分别为Ⅰ，Ⅱ，Ⅲ，Ⅳ，Ⅴ，数字越小的围岩性质越好；

Ⅰ类：岩石新鲜完整、构造影响轻微、节理裂隙不发育或稍发育，闭合且延伸不长，无或很少软弱结构面、断层带宽[0，0.1)米，与洞向近正交、岩体呈整体或块状砌体结构；

Ⅱ类：岩石新鲜或微风化，受构造影响一般，节理裂隙稍发育或发育，有少量软弱结构面、层间结合差；断层破碎带宽[0.1，0.5)米、与洞向斜交或正交、岩体呈块状砌体或层状砌体结构；

Ⅲ类：岩石微风化或弱风化，受地质构造影响裂隙发育、部分张开充泥，软弱结构面分布较多、断层破碎带[0.5，1)米，与洞线斜交或平行、岩石呈碎石状镶嵌结构；

Ⅳ类：与III类同，断裂及软弱结构面较多，断层破碎带[1，2)，与洞平行，岩体呈碎石状镶嵌结构，局部呈碎石状压碎结构；

Ⅴ类：散体：砂层滑坡堆积及碎、卵、砾质土；根据围岩属性确认该里程段区域内的分形维数，分形维数为该区域内岩石围岩级别数；

步骤三、确定三维裂隙LS文法结构：

自然界中裂隙不规则，可由随机LS文法生成，这里探讨分形裂隙的LS文法构图算法；分形几何中的LS文法构图算法是仿照语言学中语法生成的发放来构造图形的一种算法；其核心概念是重写，它制定一个或几个初始字母和一组“生成规则”，将生成规则反复作用到初始字母和新生成的字母上，产生出整个语言。

字母表：L，R；生成规则：L→R，R→LR；初始字母：R；则有：

R→LR→RLR→LRRLR→……

该文法用字母表和符号串来表达生成对象的初始形式，然后根据一组产生式重写规则，将初始形式的每个字符一次替换为新的字符形式，一次过程反复替换重写，最后生成终极图形；

LS文法图形的生成过程是将绘图的每一时刻状态定义为当前位置矢量T与前进方向a的集合(T、a)，方向a共包括裂隙生长的倾向α和倾角β两部分。根据实际工程中裂隙的调查统计结果，一般情况下，定义岩体裂隙的三种生长方式：

P₁:F→F[+F]F[-F]F，

P₂:F→F[+F]F[-F[+F]]，

P₃:F→FF[-F+F+F]+[+F-F-F]；

以第一种生长方式为例，将文法字母表的绘图规则做些简要的说明，其分形元的绘图过程以及具体符号字母的含义如下：“F”完成第一条线段AB→“[”处，程序自动记录点B坐标、方向信息→“+”表示划线的方向开始顺时针旋转预定角度a→“F”完成第二条线段BC→

“]”处，自动读取点B信息，并且返回至点B→F完成第三条线段BD→“[”记录点D信息→“-”划线方向开始逆时针旋转a→“F”完成第四条线段DE→“]”选取点D信息，并返回点D→“F”完成第五条线段DF；

此后，结果迭代，F[+F]F[-F]F中的每个F都会被F[+F]F[-F]F所替代，那么绘制的图形就会变成F[+F]F[-F]F[+F[+F]F[-F]F]F[+F]F[-F]F[-F[+F]F[-F]F]F[+F]F[-F]F；像这样一直迭代下去，每一次替换，相当于分形裂隙的一次生长；

同理，第二种生长方式和第三种生长方式与第一种生长方式类似，其基本绘制过程是：

第二种生长方式：A→B→C回到B→D→E→F回到D；

第三种生长方式：A→B→C→D→E→F回到C→G→H→I回到C；

在编制程序代码绘制分形裂隙时，需要利用随机数函数生成随机数，该随机数函数将会给p一个新的随机的数值，这样一来，p值将会被固定，就不能采用随机的生长方式；

因为在定义生长方式的时候定义的是三种生长方式，所以在编写代码时将[0，1]划分为3个区间，然后由p值所处于的区间段选取生长方式，来完成分形裂隙的生长、绘制，多次运行程序所生成的分形裂隙形状一般不同，因为从最开始生长一次时分形裂隙的形状就有3种可能，随后又因随机的p的值不同采取不同的生长方式和顺序，所以才会模拟出不同姿态的分形裂隙；

该文法尽管可以比较逼真地模拟出树的形状，但是根据文法可以看出：在设计分形元的时候，生长角度是被赋予了一个定值，这样一来，分形元的形状就被固定了下来，如果绘制出来的分形裂隙需要调整的话，就必须回到程序里面，重新编写程序，改写生长角度，这样就显得很麻烦；

另外一方面，在选择生长方式的时候，通常是将[0，1]划分为[0，0.333]、(0.333，0.666]和(0.666，1]三个均匀的区间；此时，各种生长方式被选中的概率是相同的；它所绘制的分形裂隙树形比较随机，不同时刻绘制的分形裂隙差别较大，各种生长方式被选择的概率不受人的控制；

在在绘制分形裂隙的时候，如果能够将生长元中的生长角度参数化，直接在需要调整的时候改变参数，从而实现树形的改变，将会大大提高绘图速度，改善绘图效果；在编制程序进行绘图的时候，添加一个文本框，用于输入a的值，这样便于调整生长角度，从而实现分形裂隙绘制的参数化；

LS文法中的概率p，虽然是一个随机的变量，从表面上看，分形裂隙的绘制和p的值关系不大，但是在实际的绘图过程中，概率p却起到很重要的作用；如果能控制p在不同的区间选择不同的生长方式，就可以对分形裂隙中各种生长元的比例起到人为干预作用，从而使绘制出来的分形裂隙能够更加贴近人的意愿；

为了控制树形，再引入参数U和V，这样可以将[0,1]划分为[0,U]、(U,V]和(V,1]三个区间，取不同的U值和V值，就可以控制小区间的大小，从而改变三种生长方式被选中的机率；

步骤四、重取样生成的三维模拟裂隙：

在隧道工程区不同里程段生成的三维裂隙在空间上未加约束限定条件，部分裂隙生长到地表以上，与实际情况不符，于是需要对生成的三维空间裂隙进行重采样；本发明中重取样生成的三维裂隙的解决方案是对整体区域工程区地质体有限元化，然后判断生成的分形裂隙节点是否在细化的有限元单元内部，若在内部，则保留该分形单元节点，若不在内部则去除该裂隙节点；

区域工程区地质体有限元的思路是先对隧道工程区DEM表面栅格化，地表共有点p₀，p₁，p₂，p₃，……，p_n，坐标分别为：(x₀,y₀,z₀)，(x₁,y₁,z₁)，(x₂,y₂,z₂)，(x₃,y₃,z₃)，…，(x_n,y_n,z_n)。区域工程区地表点密集，x方向最大最小值为x_max，x_min，y方向最大最小值为y_max，y_min。地表栅格边长d计算公式为：

X方向栅格个数为cx，Y方向栅格个数为cy：

cx＝(x_max-x_min)/d+1

cy＝(y_max-y_min)/d+1

第i行第j列的栅格结点坐标为(x_i，j，y_i，j，z_i，j)＝(x_min+d*i，y_min+d*j,z)，D为所有地表点到待插样点的最短距离。z值即为该处栅格点的高程：

计算每个栅格中心点坐标，第i行第j列所在栅格的中心点坐标为

第i行第j列栅格向下共有n_i，j个有限元单元：区域地质体共划分为N个有限元单元：

所有生成的三维裂隙节点坐标为(x^LS,y^LS,z^LS)，第i行第j列第k个有限元单元的中心点坐标为(0≤i≤cx,0≤j≤cy,0≤k≤n_i,j)，如果存在条件

则该三维裂隙节点位于有限元单元内部，保留该裂隙结点，否则删除该结点。

附图说明

图1为本发明的方法流程图；

图2为分形中LS文法第一种生长方式；

图3为分形中LS文法第二种生长方式；

图4为分形中LS文法第三中生长方式；

图5为生成的三维裂隙网络真实感模拟图。

具体实施方式

实施例1

本实施例提供一种基于分形算法的区域地质体岩体裂隙模拟方法，方法流程如图1所示，首先根据隧道围岩的属性信息将隧道沿线区域地质体分为不同的区域；然后对不同的区域分别做研究，根据不同区域的围岩属性信息计算此区域分形深度；之后根据实地统计信息确定三维裂隙的LS文法结构；最后重取样生成的三维模拟裂隙，可生成位于区域地质体内部真实感模拟的三维岩体裂隙。

以跃龙门隧道路线为例，包含以下步骤：

步骤一、根据隧道围岩属性划分不同区域，各结点坐标如下：

该隧道共有24个结点，为方便数值计算，此处围岩分界点处的里程数据减去隧道起始点的里程数据，准换为距离数据，共有99个围岩分界点；

围岩分界点差值里程：

根据等式：

x_i＝x_j+(x_j+1-x_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)

y_i＝y_j+(y_j+1-y_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)

z_i＝z_j+(z_j+1-z_j)·(l_i-L_j)/(L_j+1-L_j)(L_j<l_i<L_j+1)(0<i<m,0<j<n)

计算各围岩分界点的三维空间坐标为：

ID，X，Y，Z

1，9883.78，8109.41，981.5

2，10007.1，7957.07，982.737

3，10227.3，7685.04，984.946

4，10397.1，7475.19，986.899

5，10453.8，7405.24，988.541

6，10673.9，7133.21，994.928

7，10891.9，6822.03，1001.86

8，11247.7，6208.21，1014.82

9，11272，6164.54，1015.73

10，11296.4，6120.87，1016.64

11，11320.7，6077.2，1017.55

12，11534.9，5692.9，1025.58

13，11559.3，5649.23，1026.49

14，11598.2，5579.36，1027.95

15，11622.6，5535.69，1028.86

16，11744.3，5317.34，1033.43

17，11768.6，5273.67，1034.34

18，11793，5230，1035.25

19，11817.3，5186.33，1036.16

20，11866，5098.99，1037.99

21，11890.3，5055.32，1038.9

22，12046.1，4775.83，1044.74

23，12070.5，4732.16，1045.65

24，12259.4，4414.23，1052.4

25，12414，4181.01，1057.51

26，12463.2，4117.95，1058.97

27，12771.8，3777.57，1067.36

28，12806.5，3741.53，1068.28

29，12908.9，3646.56，1070.83

30，12948.1，3615.48，1071.74

31，13010.7，3565.75，1073.2

32，13049.9，3534.67，1074.11

33，13089，3503.59，1075.03

34，13128.2，3472.51，1075.94

35，13167.4，3441.44，1076.85

36，13206.5，3410.36，1077.76

37，13246.5，3380.44，1078.68

38，13315.4，3339.74，1080.14

39，13358.4，3314.31，1081.05

40，13401.5，3288.88，1081.96

41，13444.5，3263.44，1082.87

42，13899，3034.33，1092.18

43，14323.8，2887.59，1100.39

44，14401.5，2868.67，1101.85

45，14577.2，2829.6，1105.13

46，14865.1，2795.06，1110.42

47，15023.9，2776.01，1113.34

48，15093.9，2775.2，1114.62

49，15273.9，2775.82，1117.9

50，15433.8，2776.54，1120.82

51，15652.6，2799.7，1124.84

52，15870.1，2832，1128.85

53，15967.5，2854.57，1130.68

54，16335.3，2949.76，1137.61

55，16661.2，3046.68，1143.81

56，16709.1，3060.93，1144.73

57，16757，3075.18，1145.64

58，16804.9，3089.43，1146.55

59，16852.9，3103.68，1147.46

60，17226.6，3214.85，1154.58

61，17725，3363.07，1164.07

62，17792.1，3383.02，1165.34

63，18022.2，3451.43，1169.72

64，18280.9，3528.39，1174.65

65，18491.8，3591.1，1178.66

66，18549.3，3608.2，1179.76

67，18731.4，3662.36，1183.23

68，18779.3，3676.61，1184.14

69，18884.8，3707.96，1186.14

70，19172.3，3793.47，1191.62

71，19316.1，3836.23，1194.36

72，19651.5，3935.99，1200.74

73，19795.3，3978.75，1203.48

74，20226.6，4107.02，1211.69

75，20782.5，4272.34，1222.27

76，20830.4，4286.59，1223.18

77，21079.6，4360.7，1227.93

78，21233，4406.31，1230.85

79，21980.5，4628.64，1245.08

80，22076.4，4657.15，1246.9

81，22488.5，4779.72，1254.75

82，22536.4，4793.98，1255.66

83，23140.3，4973.57，1267.16

84，23303.2，5022.03，1270.26

85，23945.3，5213.02，1282.48

86，24022，5235.83，1283.94

87，24137，5270.03，1286.13

88，24213.7，5292.84，1287.59

89，24740.9，5449.62，1297.63

90，24836.7，5478.13，1299.45

91，25085.3，5554.34，1304.2

92，25161.6，5578.23，1305.66

93，25314.3，5626.01，1308.57

94，25446.5，5671.76，1311.13

95，25905，5869.99，1320.25

96，25949.6，5892.58，1321.16

97，26065.6，5951.33，1323.54

98，26110.2，5973.93，1324.45

99，26476.3，6198.59，1327.11

步骤二、根据不同区域的围岩属性信息计算此区域分形深度：

隧道围岩属性根据破碎程度及岩性分为5级，分别为Ⅰ，Ⅱ，Ⅲ，Ⅳ，Ⅴ，数字越小的围岩性质越好；

Ⅰ类：岩石新鲜完整、构造影响轻微、节理裂隙不发育或稍发育，闭合且延伸不长，无或很少软弱结构面、断层带宽[0，0.1)米，与洞向近正交、岩体呈整体或块状砌体结构；

Ⅱ类：岩石新鲜或微风化，受构造影响一般，节理裂隙稍发育或发育，有少量软弱结构面、层间结合差；断层破碎带宽[0.1，0.5)米、与洞向斜交或正交、岩体呈块状砌体或层状砌体结构；

Ⅲ类：岩石微风化或弱风化，受地质构造影响裂隙发育、部分张开充泥，软弱结构面分布较多、断层破碎带[0.5，1)米，与洞线斜交或平行、岩石呈碎石状镶嵌结构；

Ⅳ类：与Ⅲ类同，断裂及软弱结构面较多，断层破碎带[1，2)，与洞平行，岩体呈碎石状镶嵌结构，局部呈碎石状压碎结构；

Ⅴ类：散体：砂层滑坡堆积及碎、卵、砾质土；根据围岩属性确认该里程段区域内的分形维数，分形维数为该区域内岩石围岩级别数；

根据不同区域的围岩属性信息计算此区域分形深度，分形维度正相关于围岩级别数字：

步骤三、确定三维裂隙LS文法结构：

将绘图的每一时刻状态定义为当前位置矢量T与前进方向a的集合(T、a)，方向a共包括裂隙生长的倾向α和倾角β两部分。根据实际工程中裂隙的调查统计结果，一般情况下，定义岩体裂隙的三种生长方式：

P₁:F→F[+F]F[-F]F，

P₂:F→F[+F]F[-F[+F]]，

P₃:F→FF[-F+F+F]+[+F-F-F]；

根据跃龙门隧道工程区的裂隙统计数据，图2、图3与图4为区域工程区地质体的裂隙微小单元结构；

生成的裂隙结点三维空间坐标为：

ID，X，Y，Z

1，10221.7，7691.93，984.89

2，10221.7，7691.93，994.89

3，10218.4，7693.84，1004.13

4，10221.7，7691.93，994.89

5，10221.7，7691.93，1004.89

6，10225，7693.84，1014.13

7，10221.7，7691.93，1004.89

8，10221.7，7691.93，1014.89

9，10218.4，7693.84，1024.13

10，10214.8，7699.97，1031.2

11，10218.4，7693.84，1024.13

12，10215.1，7695.76，1033.37

13，10215.1，7695.76，1043.37

14，10211.8，7697.67，1052.61

15，10221.7，7691.93，1014.89

16，10221.7，7691.93，1024.89

17，10218.4，7693.84，1034.13

18，10221.7，7691.93，1024.89

19，10221.7，7691.93，1034.89

20，10225，7693.84，1044.13

21，10221.7，7691.93，1034.89

22，10221.7，7691.93，1044.89

23，10225，7693.84，1054.13

24，10225，7693.84，1064.13

25，10225，7693.84，1054.13

26，10228.3，7695.76，1063.37

27，10231.9，7701.88，1070.44

28，10235.2，7703.79，1079.68

29，10221.7，7691.93，1044.89

30，10221.7，7691.93，1054.89

31，10218.4，7693.84，1064.13

32，10221.7，7691.93，1054.89

33，10221.7，7691.93，1064.89

34，10225，7693.84，1074.13

35，10221.7，7691.93，1064.89

36，10221.7，7691.93，1074.89

37，10218.4，7693.84，1084.13

38，10214.8，7699.97，1091.2

39，10218.4，7693.84，1084.13

40，10215.1，7695.76，1093.37

41，10215.1，7695.76，1103.37

42，10215.1，7695.76，1093.37

43，10211.8，7697.67，1102.61

44，10208.4，7699.58，1111.85

45，10204.9，7705.71，1118.92

46，10208.4，7699.58，1111.85

47，10205.1，7701.5，1121.08

48，10205.1，7701.5，1131.08

49，10201.8，7703.41，1140.32

50，10205.1，7701.5，1121.08

51，10205.1，7701.5，1131.08

52，10205.1，7701.5，1141.08

53，10208.4，7703.41，1150.32

54，10208.4，7703.41，1160.32

55，10205.1，7705.32，1169.56

56，10205.1，7701.5，1141.08

57，10201.6，7707.62，1148.16

58，10198.3，7709.53，1157.39

59，10198.3，7709.53，1167.39

60，10205.1，7705.32，1169.56

61，10201.6，7711.45，1176.63

62，10201.6，7720.69，1180.46

63，10201.6，7711.45，1176.63

64，10198.1，7717.57，1183.7

65，10198.1，7717.57，1183.7

66，10194.5，7723.7，1190.78

67，10201.8，7703.41，1140.32

68，10201.8，7712.65，1144.15

69，10206.8，7721.31，1144.15

70，10201.8，7712.65，1144.15

71，10201.8，7721.89，1147.98

72，10198.3，7728.01，1155.05

73，10201.8，7721.89，1147.98

74，10201.8，7731.13，1151.8

75，10198.3，7737.25，1158.87

76，10198.3，7746.49，1162.7

77，10198.3，7737.25，1158.87

78，10194.7，7743.38，1165.95

79，10191.4，7745.29，1175.18

80，10187.9，7751.41，1182.26

81，10194.7，7743.38，1165.95

82，10191.4，7745.29，1175.18

83，10188.1，7747.2，1184.42

84，10188.1，7747.2，1184.42

85，10188.1，7756.44，1188.25

86，10221.7，7691.93，1074.89

87，10221.7，7691.93，1084.89

88，10218.4，7693.84，1094.13

89，10221.7，7691.93，1084.89

90，10221.7，7691.93，1094.89

91，10225，7693.84，1104.13

92，10221.7，7691.93，1094.89

93，10221.7，7691.93，1104.89

94，10218.4，7693.84，1114.13

95，10214.8，7699.97，1121.2

96，10218.4，7693.84，1114.13

97，10215.1，7695.76，1123.37

98，10215.1，7695.76，1133.37

99，10215.1，7695.76，1123.37

100，10211.8，7697.67，1132.61

101，10221.7，7691.93，1104.89

102，10221.7，7691.93，1114.89

103，10221.7，7691.93，1124.89

104，10225，7693.84，1134.13

105，10225，7693.84，1144.13

106，10221.7，7695.76，1153.37

107，10221.7，7691.93，1124.89

108，10218.2，7698.05，1131.96

109，10214.8，7699.97，1141.2

110，10214.8，7699.97，1151.2

111，10221.7，7695.76，1153.37

112，10221.7，7695.76，1163.37

113，10221.7，7695.76，1173.37

114，10225，7697.67，1182.61

115，10221.7，7695.76，1173.37

116，10218.2，7701.88，1180.44

117，10214.8，7703.79，1189.68

118，10214.8，7699.97，1151.2

119，10214.8，7699.97，1161.2

120，10211.5，7701.88，1170.44

121，10214.8，7699.97，1161.2

122，10214.8，7699.97，1171.2

123，10218.2，7701.88，1180.44

124，10214.8，7699.97，1171.2

125，10214.8，7699.97，1181.2

126，10211.5，7701.88，1190.44

127，10214.8，7699.97，1181.2

128，10214.8，7699.97，1191.2

129，10334，8051.17，1188.79

130，10334，8051.17，1178.79

131，10334，8054.99，1169.55

132，10334，8051.17，1188.79

133，10337.5，8057.29，1181.72

134，10345.5，8061.91，1177.89

135，10337.5，8057.29，1181.72

136，10341，8063.41，1174.65

137，10341，8067.24，1165.41

138，10341，8063.41，1174.65

139，10344.6，8069.54，1167.58

140，10334，8054.99，1169.55

141，10334，8058.82，1160.32

142，10334，8062.65，1151.08

143，10334，8062.65，1141.08

144，10334，8066.47，1131.84

145，10337.5，8072.6，1124.77

146，10334，8062.65，1151.08

147，10342，8067.27，1147.25

148，10345.5，8073.39，1140.18

149，10345.5，8077.22，1130.94

150，10337.5，8072.6，1124.77

151，10337.5，8076.43，1115.53

152，10337.5，8080.25，1106.29

153，10337.5，8080.25，1096.29

154，10337.5，8084.08，1087.05

155，10341，8090.2，1079.98

156，10337.5，8080.25，1106.29

157，10345.5，8084.87，1102.46

158，10349，8091，1095.39

159，10349，8094.82，1086.15

……

……

……

步骤四、重取样生成的三维模拟裂隙：

在隧道工程区不同里程段生成的三维裂隙在空间上未加约束限定条件，部分裂隙生长到地表以上，与实际情况不符，于是需要对生成的三维空间裂隙进行重采样；本发明中重取样生成的三维裂隙的解决方案是对整体区域工程区地质体有限元化，然后判断生成的分形裂隙节点是否在细化的有限元单元内部，若在内部，则保留该分形单元节点，若不在内部则去除该裂隙节点；

区域工程区地质体有限元的思路是先对隧道工程区DEM表面栅格化，地表共有点p₀，p₁，p₂，p₃，……，p_n，坐标分别为：(x₀,y₀,z₀)，(x₁,y₁,z₁)，(x₂,y₂,z₂)，(x₃,y₃,z₃)，…，(x_n,y_n,z_n)。区域工程区地表点密集，x方向最大最小值为x_max，x_min，y方向最大最小值为y_max，y_min。地表栅格边长d计算公式为：

X方向栅格个数为cx，Y方向栅格个数为cy：

cx＝(x_max-x_min)/d+1

cy＝(y_max-y_min)/d+1

第i行第j列的栅格结点坐标为(x_i，j，y_i，j，z_i，j)＝(x_min+d*i，y_min+d*j,z)，D为所有地表点到待插样点的最短距离。z值即为该处栅格点的高程：

计算每个栅格中心点坐标，第i行第j列所在栅格的中心点坐标为

第i行第j列栅格向下共有n_i,j个有限元单元：区域地质体共划分为N个有限元单元：

所有生成的三维裂隙节点坐标为(x^LS,y^LS,z^LS)，第i行第j列第k个有限元单元的中心点坐标为(0≤i≤cx,0≤j≤cy,0≤k≤n_i,j)，如果存在条件

则该三维裂隙节点位于有限元单元内部，保留该裂隙结点，否则删除该结点，最后生成的裂隙网络如图4所示。

以上实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。