基于等高线的复杂地形三维离散单元法数值模型建模方法

摘要

本发明提出的基于等高线的复杂地形三维离散单元法数值模型建模方法，属于地质工程及采矿工程领域中数值模拟研究领域。该方法首先获取建模区域的数字地形图，在区域边界和高程点稀疏的区域加密高程点并求取每个加密高程点的高程值，提取每个高程点的三维空间信息并保存；然后对数字地形图进行地表曲面样条插值，并绘制建模区域地表三维图形；将曲面样条插值后的所有高程点的三维空间信息输出到最终地形数据文件中，通过三维离散单元法软件建立三维数值模型；通过与地表三维图形对比验证，得到建模区域最终的三维离散单元法数值模型。本发明可对任一复杂地形区域进行高精度的数值建模，所建模型地表曲面与实际地形表面偏差小，有较高的应用价值。

说明书

技术领域

本发明属于地质工程及采矿工程领域中的数值模拟研究领域，具体地涉及一种基于等高线的复杂地形三维离散单元法(3DEC)数值模型建模方法。

背景技术

地形地貌是普通地图上需要表示的重要基本要素，通常也将其称为地势或地形，特指地表的高低起伏形态。地形在地图上的表示方式有多种，如写景法、晕渲法、分层设色法和等高线法。其中等高线法最为常用，且表示精度最高，它的基本原理是：将地面高程相等的点连成平滑的曲线并将其投影到平面上。等高线在实际工程建设中有多种用途，如两点之间坡度的计算，工程开挖土方量的计算，水利工程汇水面积的计算，地表任一点海拔的确定等。

在岩土工程和采矿工程的数值模拟研究中，三维有限差分程序FLAC3D(FastLagrangian Analysis of Continua)和三维离散单元法程序3DEC(3Dimension DistinctElement Code)是最主要的两种软件，在进行数值模拟之前，首先应建立与模拟区域地形及地下岩层相一致的数值计算模型。当前，由于建模技术方法的限制，很多研究者在数值模拟计算中较少考虑实际地形起伏状况对数值模拟结果的影响，当然，对于平原地区，这种近似的建模方法是可行的，但在丘陵区或山区，由于地形起伏较大，如果不考虑地形变化的影响，模型过于简化，则数值模拟结果与现场实测结果将会存在较大误差，严重制约数值模拟计算的精度及其结果的可靠性。

针对这一问题，尽管有的研究者对复杂地形的数值建模进行了一定程度的研究，如尹士献博士在公开发表的论文《复杂地形三维数值建模研究》中提出：对于简单地形要采用空间曲面方程进行描述，对于复杂地形必须首先拟合出空间曲面方程，然后再利用数值软件中自定义函数功能建立三维数值模型。由于并不是所有的地形均可用曲面方程的形式加以描述，因此尹博士提出的方法，对于地形并不十分复杂的地区是较为适用的，但对于地形复杂的山区，这种方法便存在一定的局限性，建模精度将会大大降低，难以满足较高精度数值模拟计算的要求，另外，该方法前期数据处理工作量大，在实际操作时过程也较为复杂。

发明内容

本发明的目的是为克服已有技术的不足之处，提出一种基于等高线的复杂地形三维离散单元法数值模型建模方法。本发明所给出的建模方法可对任一复杂地形区域进行高精度的数值建模，智能化程度高，所建数值模型地表曲面与实际地形表面偏差很小，有较高的实际应用价值。

本发明提出的基于等高线的复杂地形三维离散单元法数值模型建模方法，其特征在于，包括以下步骤：

1)获取建模区域的数字地形图；

确定研究区并获取研究区的数字地形图；对研究区数字地形图进行裁剪得到建模区域的数字地形图，所述建模区域的数字地形图上包括若干高程点；

2)在步骤1)得到的建模区域数字地形图上的，在区域边界和高程点稀疏的区域加密高程点，然后利用解析法求取每个加密高程点的高程值；

3)利用步骤2)完成加密高程点后的建模区域数字地形图，提取每个高程点的三维空间信息：包括每个高程点的(x,y)坐标和高程值，将所有高程点的三维空间信息保存为TXT或XLSX形式的建模区域地形数据文件；

4)对步骤3)完成后的建模区域数字地形图进行地表曲面样条插值，得到经曲面样条插值后的所有高程点的三维空间信息，并利用MATLAB绘制建模区域地表三维图形；

5)利用建模区域的最终地形数据文件，建立该区域的三维数值模型；

将经过步骤4)曲面样条插值后的所有高程点的三维空间信息，通过MATLAB语言程序模块按照三维离散单元法软件可识别的命令格式输出到一个TXT文件中，将该TXT文件记为建模区域的最终地形数据文件；打开三维离散单元法软件，导入该最终地形数据文件并执行，建立该建模区域的三维数值模型；

6)通过对比验证，得到建模区域最终的三维离散单元法数值模型；

将步骤5)建立的三维数值模型和步骤4)用MATLAB绘制的地表三维图形进行对比并判定：

如果步骤5)得到的三维数值模型的上表面与MATLAB所绘制的地表三维图形不一致，则说明步骤5)中输出的建模区域的最终地形数据文件有误，重新返回步骤5)，对MATLAB语言程序模块进行修改，然后生成新的三维数值模型；否则，步骤5)得到的三维数值模型即为该建模区域的最终三维离散单元法数值模型。

本发明的特点及有益效果在于：

本发明提出的提出一种基于等高线的复杂地形三维离散单元法数值模型建模方法，该方法的优点主要包括：一，充分利用了等高线所表达的地形信息数据，并可根据建模精度的需要确定高程点的加密间距；二，获取地形信息后，采用MATLAB编程可实现地表三维图形的自动绘制，并同时输出3DEC可识别的命令文件；三，将命令文件导入3DEC软件中，可自动实现复杂地形区域的三维数值建模，从而大大提高了建模的智能化程度，且通过将3DEC数值模型与MATLAB所绘三维图形进行对比，可检验建模结果的可靠性。

利用本发明所描述的方法建立模型，其前期数据处理工作量小、智能化程度高，建模速度快，且能够建立任意复杂地形的三维数值模型，所建数值模型地表曲面与实际地形表面偏差很小，并且经过微小改动还可推广应用到FLAC3D数值建模中。在岩土工程和采矿工程领域数值模拟研究中，采用本发明所描述的方法进行数值建模，不但可以大大节省研究人员的建模工作量，而且不用再担心实际工程区域地形的复杂情况，可以完全按照实际地形进行建模，且建模精度较高，如果能够再合理的选择地层相应的岩性参数，那么数值模拟计算结果与工程实际之间的差距将会大大减小，从而提高利用数值模拟解决实际工程问题的可靠性，推动数值模拟技术在工程中的广泛应用。

附图说明

图1是本发明方法的流程图。

图2是本发明实施例中建模区域数字地形图的示意图。

图3是本发明实施例中所要提取的高程点及加密高程点示意图。

图4本发明实施例中样条插值后绘制的地表三维图形示意图。

图5是本发明实施例的3DEC数值模型示意图。

具体实施方式

本发明提出的基于等高线的复杂地形三维离散单元法数值模型建模方法，下面结合附图和具体实施例对本发明做进一步详细说明。

本发明提出的基于等高线的复杂地形三维离散单元法数值模型建模方法，整体流程如图1所示，包括以下步骤：

1)获取建模区域的数字地形图；

确定研究区并获取研究区的数字地形图；研究区数字地形图是获得地表点三维坐标数据的基础文件，；假定要进行某区域的数值模型创建，那么首先要收集该区域的数字地形图，数字地形图的文件名一般是以.DWG为后缀；数值地形图通常由甲方提供或委托测绘单位进行现场测绘。

通常收集到的研究区数字地形图涉及的范围会比所要建模的区域大很多，为了提取数据的方便，需要将收集到的数字地形图进行裁剪，或删掉建模区域以外的所有地物和地形信息。根据研究区中实际建模区域的大小对研究区数字地形图进行裁剪得到建模区域的数字地形图(为了提取信息不受干扰，通常需要删除数字地形图上地表点三维坐标以外的地物信息)。

本实施例建模区域数字地形图的示意图如图2所示，图中的点为高程点，数字代表的该点的高程值。

2)在步骤1)得到的建模区域数字地形图上加密高程点；加密完成后，利用解析法获得每个加密高程点的高程值；

由于建模区域四周的高程点可能比较稀疏，为提高所建模型边缘区域的精度，需在建模区域的边界周围加密一定数量的高程点(高程点加密的数量没有具体的限制，通常加密的高程点越多所建模型精度越高)，本实施例的所要提取的高程点及加密高程点示意图如图3所示；图3中，旁边附有数值的点为原数字地形图上的加密高程点，矩形区域边缘的高程点(没有附上数值的高程点)即为加密高程点，这些点的坐标需要在数字地形图(在原地形图或建模区域地形图都可以)上进行提取，而高程值则需要通过解析法计算获得。另外，如果地形图上还存在高程点稀疏的区域，或是为了更好的表达一些特殊的地貌特征，如山谷、山脊等，可在数字地形图中的相关区域加密高程点，高程点的加密方法采用解析法。

3)利用步骤2)完成加密高程点后的建模区域数字地形图进行三维空间信息提取，并保存为建模区域的地形数据文件；

当建模区域高程点加密完成之后，即可采用AutoCAD软件中的数据提取工具，获得建模区域每个高程点的三维空间信息，包括每个高程点的点位坐标(x,y)和高程值H；提取后将所有高程点的三维空间信息保存为TXT或XLSX形式的建模区域地形数据文件。

利用AutoCAD软件提取高程点，操作流程如下：打开软件，在工具栏里点击“数据提取”菜单，在弹出的提示窗里选择“创建新数据提取”，然后点击“下一步”，这时会弹出对话框，提示输入文件名，这个名字可任意输入，保存的位置也可任意选择，此文件只是一个记录提取规则的形式文件，并不是最终需要的数据文件；文件名输入后，点击“保存”按钮，系统会弹出对话框，提示选择提取数据源，这里有两个选择，即：选择“包括当前图形”还是选择“在当前图形中选择对象”，因为在第一步已将其他不相关的信息删除了，所以在这里直接选择后者；再点击“下一步”，在弹出的对话框中选择高程点所在的图层，本实施例中采用南方CASS制图软件所制作的地形图，其高程点所在图层为“GCD”，直接选择即可，其他软件制作的地形图，如果高程点没有Z属性，只有文字，那么在提取时只提取文字的值和文字的x和y坐标，即为高程点信息；接着点击“下一步”，此时可以看到所选对象包含的所有特性，选择“几何图形”属性里的X、Y、Z值即可，继续下一步，便可观察到所提取的数据，将提取出的数据输出到EXCEL，也可输出到TXT文件，为便于后续的数据处理，通常将其输出到EXCEL文件。

步骤3)的作用是生成MATLAB软件可以直接编程读取的数据文件，为后续三维绘图和数值建模做准备。步骤3)是进行三维绘图的基础数据，通常情况下通过地形图获取的高程点及加密高程点数量还是太少，还不足以绘制精度较高三维地形，还需要进行样条差值，因此步骤3)是步骤4)的基础。

4)对步骤3)完成后的建模区域数字地形图进行地表曲面样条插值，得到经曲面样条插值后的所有高程点的三维空间信息，并利用MATLAB绘制建模区域地表三维图形；

经过步骤2)加密高程点之后，建模区域数字地形图已经具备了建模区地表点的基础三维坐标数据，但对于绘制高精度的三维地形表面而言，经过前述步骤提取的高程点和加密的高程点信息仍满足不了需要，点的数量还很不足，需要采用曲面样条插值的方法对高程点继续加密，在进行曲面样条插值时可以人为的设定插值的x方向距离和y方向距离，以控制模型的精度。插值后采用MATLAB调取数据编程，绘制加密后的地表三维图形，本实施例经过样条插值后绘制的三维地表三维图形如图4所示。图4中，坐标轴分别是：x轴，y轴和z轴，色块是由MATLAB程序自动生成的，以使三维图看起来更逼真。

曲面样条插值方法具体步骤如下：

加密前的高程点坐标和高程值可采用式(1)表示:

式中：x_i为高程点i的x坐标，y_i为高程点i的y坐标，H_i为高程点i的高程值。

定义公式(1)为二元单值列表函数，进而对该函数进行拟合，其二元样条函数表达式如(2)所示：

式中：r²＝(x-x_i)²+(y-y_i)²，c₁、c₂、…、c_3+n为待求系数，ε为调节系数，一般采用经验值，对于平坦地区ε介于0.01～1之间，对于奇异曲面ε介于10^-6～10^-5之间；

待求系数则可以由下式确定：

式中：r_ji²＝(x_j-x_i)²+(y_j-y_i)²；h_j为第i个节点的加权系数，通常情况下在一般的插值计算时，h_j的取值都可为0，以使所拟合的曲面与给定原始高程点数据相吻合。如果将式(3)用矩阵的形式表示，则表达式如式(4)所示：

A_m×mc_m×1＝H_m×1(4)

式(4)中的系数矩阵c_m×1是由节点坐标值和加权系数所组成的对称阵，如果A_m×m不是奇异矩阵，则可求解方程，进而得出系数矩阵，表达式如式(5)所示：

系数矩阵求出之后，式(2)则同时被确定，如果将拟合区域(即建模区域)划分为m个四边形网格，那么每个网格节点所对应的高程可用函数(6)表示：

由于任一网格节点x_k处的高程函数对坐标x和y一阶偏导数函数如式(7)所示：

将公式(7)与四边形网格的节点坐标相结合，可求得拟合面网格节点所对应的空间坐标，从而得到经过曲面样条插值加密之后的高程点三维坐标。

5)利用建模区域的最终数据地形文件，建立该区域的三维数值模型；

用MATLAB语言编写程序模块，将经过曲面样条插值后的建模区域数字地形图上所有高程点三维空间信息，按照3DEC软件可识别的命令格式输出到一个TXT文件中，将该TXT文件记为建模区域的最终地形数据文件。打开数值模拟软件3DEC，导入输出最终数据地形文件的TXT文件，执行该文件，即可建立相应的建模区域的三维数值模型，

3DEC数值建模的关键是，将经过插值之后的每一个四边形网格点的空间坐标(x，y，H)，按照先纵向后横向或先横向后纵向的方式以3DEC软件相应绘图命令prism或face所要求的格式进行输出，输出的格式通常为TXT格式；

本实施以3DEC内置建模命令prism为例，说明文件的输出格式。

3DEC poly prism命令流格式如下

操作时，打开3DEC软件，在文件中选择open new item，选择输出的TXT数据文件，执行之后，即可生成如图5所示的数值模型。

经过上述各步骤，图5所建立的数值模型即为根据图2中的地形数据生成的，采用三维图形进行表达，三维坐标系的方向如图中左下角的坐标轴所示，如果将该三维图形的上表面向下垂直投影，则投影的图形应与图4相似，相似度越高，说明建模的精度越高。

6)通过对比验证，得到建模区域最终的三维离散单元法(3DEC)数值模型；

将步骤5)建立的三维数值模型和步骤4)用MATLAB绘制的地表三维图形进行对比，也可结合图2进行对比，检验三维数值模型的正确性：

如果步骤5)得到的三维数值模型的上表面与MATLAB所绘制的地表三维图形不一致，则说明在步骤5)中输出的建模区域的最终地形数据文件不正确，此时应检查MATLAB语言程序模块是否存在错误，如果有错误则需要修改程序，重新返回步骤5)，然后重新生成新的三维数值模型；否则，步骤5)得到的三维数值模型即为该建模区域的最终3DEC数值模型，本实施例所建立的最终3DEC数值模型如图5所示。

数值模型的创建是进行数值模拟工作的前提，当模型建好之后，可根据实际工程情况对模型进行进一步的分层建模和参数赋值以满足不同工程的需要。