张 青，胡 成

(中国地质大学(武汉)环境学院，湖北 武汉 4300741)

地下热流及水流数值模拟程序Tough2发布于1991年，经过20多年的开发，已发展为现今的新流体模块，应用范围极其广阔，如CO2地质封存、地下水流模拟、环境影响评价、地热储藏工程等。虽然Tough2可以模拟各类情况，但是推广上却较为困难；需要用户具备一定的编程基础[1- 3]。在模拟过程的初期，模型构建和网格剖分工作量巨大，后人开发出相应的网格生成工具，如MeshMaker、PetraSim、MView等，但是这些网格生成工具并不具有地理信息的概念，也缺乏对网格图形强大的处理能力。

ArcGIS是具有强大空间分析与处理功能的可视化软件，可以对具有空间信息的多边形进行编辑，其平台提供的二次开发Python语言是一种跨平台的开源脚本语言，具有简单易学、开发效率高、伸缩性和可移植性高、可嵌入和稳定成熟等优势[4- 5]。自ArcGIS 9.0引入Python作为脚本语言发展至今，Python已成为一种运用到ArcGIS中实现其数据管理与分析等多种功能自动化语言，Python可以创建在独立脚本环境、工具箱或模型中运行的脚本。同时，Python可以轻松调用ArcGIS中各种工具，或访问其他人员开发的Python模块[6- 7]，从而实现程序化处理复杂数据，以减少工作量和提高工作效率。本文根据数值模拟项目实际需要，在ArcGIS平台通过Python实现Tough2多边形网格的生成，为Tough2的进一步模拟提供便捷的基础。

1 方法与原理

1.1 Tough2网格文件规则

Tough2模拟器需要的网格文件包括两部分，ELEME单元文件和CONNE连接文件，前者反映每个网格单元的基本信息，如单元标识、XYZ坐标、参数分区标识等；后者反映每个单元网格之间的关系，如交界面的面积、网格中心到交界面的距离、重力方向与网格中心连线的余弦值、热交换因子等，具体参数见表1和表2。其实，大部分参数都是网格固有的属性，但有些参数是需要根据实际情况进行赋值的，如参数分区标识和热交换因子。实际上Tough2支持任意形状的网格剖分，无论是三角性、矩形还是任意多边形都是可以的。

1.2 泰森多边形生成原理

在一定范围内生成离散点，利用离散点构建三角形网，并作出每个三角形各边的垂直平分线，每个离散点被各边垂线形成的一个多边形包围，这个多边形即为泰森多边形[8- 9]。泰森多边形网格非常符合Tough2模拟器的特点，其不规则三角网的建立是根据Delaunay准则生成，在网格剖分上更具零活性，ArcGIS平台中自带泰森多边形生成工具，利用Python语言可轻松调用该工具，并对生成的多边形进行信息提取和图形分析，如图1所示。

表1 ELEME单元文件参数说明

表2 CONNE连接文件参数说明

图1 ArcGIS随机点生成泰森多边形

1.3 网格文件信息获取

生成泰森多边形后可利用ArcGIS提供的内置函数提取网格文件所需要的坐标、体积、面积等参数，并可以利用Python基本函数计算距离、余弦值等参数，再根据实际情况赋予属性参数，最后可将提取和计算出的参数写入网格文件中。

2 技术流程

ELEME单元文件和CONNE连接文件正是Tough2所需的MESH网格文件，先有ELEME单元文件，再有CONNE连接文件；即先有网格单元，再有网格关系。本次按先获取ELEME文件再获取CONNE文件进行研究，技术流程如图2所示。

图2 技术流程图

2.1 ELEME文件参数获取

(1)随机点的生成

在网格边界内生成一些随机点，通过随机点划分出符合Delaunay准则的不规则三角网(TIN)，三角形各边的垂直平分线即可形成泰森多边形的边。在ArcGIS中提供了2种方法用于生成随机点，一种是使用创建随机点的工具，可以按照工具给出的提示生成指定范围内的随机点；另一种是调用ArcGIS中CreateRandomPoints_management()函数。在该函数中只需要给定边界文件、随机点个数、最小点间距等参数便可快速生成随机点，其函数的部分参数说明见表3。

表3 随机点函数参数说明

生成的随机点文件中，系统已经默认赋予点文件一些属性，如ID编号、X坐标、Y坐标等，但这些系统分配好的参数并不能满足Tough2网格文件生成所需要的信息，还需在随机点文件中添加材质名称这一属性。

(2)创建泰森多边形

创建泰森多边形的前提是在已知点位置的条件下，根据Delaunay准则生成的。与创建随机点类似，ArcGIS也提供了2种方法用于创建泰森多边形：一种是利用CreateThiessenPolygons工具箱；另一种是用CreateThiessenPolygons_analysis()函数，同样该函数只需要输入随机点文件和输出位置便可以生成泰森多边形，具体参数见表4。

表4 泰森多边形函数参数说明

与随机点函数不同的是，生成的泰森多边形不需要指定输出范围参数，所以，生成的泰森多边形是一个矩形，该矩形的范围是输入点输入要素的范围另加10%，在生成后还需用Clip工具剪裁至模拟区边界形状。

(3)泰森多边形网格信息提取

生成符合要求的泰森多边形后需要提取每个单元格的体积，此次研究二维网格的生成，实际上单元体积在数值上等于单元面积乘以模拟层厚度。X、Y、Z的坐标值则需要根据重力方向进行选取，默认重力方向在Z坐标轴上，竖直向下，X坐标和Y坐标等于投影坐标系减去投影带号所在范围的基数值，Z坐标等于模拟层厚度的一半。如果重力方向为Y方向，则Y坐标与Z坐标互换，Y坐标在数值上等于模拟层厚度的一半。对于平行多层网格，在空间上改变的只需是Z值。

在ArcGIS中，单元格面积和坐标是图形固有属性，上述网格信息可以通过数据访问模块(arcpy.da)中的Search Cursor类进行只读权限访问。

至此，ELEME文件的主要参数ID索引码、X坐标、Y坐标、单元体积、参数获取完成。

2.2 CONNE文件参数获取

(1)创建邻近关系表

ArcGIS中的面领域工具用于创建领近关系表，面邻域工具按照等级路径确定要在输出表中记录的邻域类型和统计数据。邻域关系、重合边和结点邻域按照从高到低的等级顺序进行重叠。一旦找到更高顺序的邻域，该工具计算和存储关系信息并跳过较低顺序关系的分析。由于在创建的泰森多边形网格(TOUGH2网格)中，领域类型不存在重叠面，需要的是重合边以获取相交面积这一参数。

在图3中，左边输入要素是一个3×3的网格单元，单元编号从101～109，右边是通过面领域工具分析后，输出的面领域分析表。

图3 面领域关系表

在图3输出的面领域关系表中，揭示了相邻网格单元间的关系，与101网格相邻的有102、104、105网格单元，102和104网格单元的领域类型是重合边，105网格单元的领域类型是结点领域。对于Tough2网格文件需要获取的是重合边，在图5中是LENGTH属性字段。输出领近关系表后通过数据访问模块(arcpy.da)中的SearchCursor访问这些信息。

(2)获取公共边到节点的距离

每两个相邻单元格都只有两个共同的端点，通过数据访问模块进行访问，比较每个网格端点的X坐标和Y坐标是否相等来获取公共端点的X坐标和Y坐标，这样每个单元格中心与两个端点间便构成了一个三角形，可通过海伦公式：

S=sqrt(p(p-a)(p-b)(p-c))p=(a+b+c)/2

和面积公式S=1/2ah求出单元格中心到公共边的距离d，即网格中心到公共边所形成等腰三角形上的高h。

实际中泰森多边形具有相邻关系的两个单元格中心到公共边的距离相等，通过两中心的坐标求出距离，其距离的一半即公共边到节点的距离。

(3)获取cosine值

在二维网格中如果取重力方向在Y坐标轴上，获取重力方向与连接线夹角的cosine值方法与计算公共边到节点的距离方式类型，以相邻网格单元的中心节点(随机点)构建直角三角形，cosine值等于两中心节点间的高差比上两中心节点的距离。在构建直角三角形前需要判断两种情况的出现，如果两节点的X坐标相等，则cosine值等于0，如果Y坐标相等，则cosine值等于1或-1。重力方向取X方向与重力方向取Y坐标轴方向相似，在此不做说明。如果重力方向在Z坐标轴上，则所有的cosine值取0。

至此，CONNE文件的主要参数相邻单元格的ID索引码、单元节点距离、cosine值的获取完成，其中，渗透性编号与网格单元的空间位置关系有关，如果两中心节点在X平面上则为1，在Y平面上则为2，在Z平面上则为3，其反映的是XYZ方向的渗透率，1、2、3只是用于识别渗透性方向的代号，可以在网格文件生成时写入。

2.3 MESH文件参数整合

在获取Tough2中的ELEME和CONNE所需的主要参数后，将这些参数按照指定的格式写入到MESH文件中，使用Python中的format函数格式化输出，例如体积的格式化输出:Volunme=format(ELEME[i][2]，′.4E′)。使用Python中的rjust可以将字符进行右对齐，如:f.write(str(D1).rjust(10))。为了将上述繁琐的步骤简单化，方便Tough2网格文件的生成，以Python工具箱的形式进行集成，其工具箱界面如图4和图5所示。

图4 MESH Maker-Polygon.pyt工具箱

图5 MESH FIle-Polygon.pyt工具箱

图6的Python工具箱用于生成指定范围内的泰森多边形网格(Tough2网格)，图7的Python工具箱用于输出MESH文件。

2.4 参数分区

由于在实际模拟项目中，研究区内各个网格单元的参数是不同的，在进行模型识别时，不可能也没有必要让每个单元都赋予不同的参数值，所以只能进行参数分区，在每个参数分区内的各单元的参数是相同的[10]。先在ArcGIS平台将研究区划分为多个子区域，每个子区域内的单元网格参数一致，在ELEME单元文件中，每一行有一个位置记录该单元的参数，用相应的代号表示，可以是英文字母也可以是数字或者组合，对于具有相同参数的单元该代号相同。在ArcGIS平台中可利用Python语言将这些落在同一个参数分区内的单元节点的参数字段统一赋予相同的代号，并对同一参数分区内的网格进行加密处理。

3 方法应用

以安徽某一石油洞库模拟为例，为了预测在库区施工及运营条件下，库区及周边地下水水位下降情况及其对地质环境可能的影响，需要在分析研究区水文地质条件的基础上，建立水文地质概念模型。根据地质时代、成因、岩性及工程性质的不同，将研究区的岩性分为三大类。其中研究区东南部主要为第四系地层，成分以粘土、砂土、砾石为主；中部主要为印支期花岗岩体，成分以钾长石、斜长石、石英、黑云母为主；西北部主要为寒武纪片岩段，成分以云母石英片岩，顶部为千枚岩为主；洞库建立在中部的花岗岩体区域。研究区参数分区如图6所示。

图6 研究区参数分区图

在ArcGIS中勾出模拟区边界的矢量文件后，在Python脚本中输入需要生成的网格数以及网格中心最小间距即可将模拟区剖分为泰森多边形网格，如图7所示。

图7 研究区生成的泰森多边形

根据落入每个参数分区内节点赋予岩性代号，并自动生成Tough2需要的MESH网格文件。本次处理的坐标系为“1954北京坐标系”，网格中心最小间距为10m，网格数量为100个，网格层厚为20m。网格文件结果如图8—9所示。

图8 生成的ELEME网格文件

图9 生成的CONNE网格文件

在图8的ELEME单元文件中，第一列代表的数值含义与表1中第一列的含义一致，即图8中的数字1代表单元标识为1的网格单元；在图9的CONNE网格文件中，第一列第一行表示相邻网格单元1和2的渗透性、距离等信息，分别与表2中的说明一一对应。通过在ArcGIS平台上生成Tough2所需的MESH网格文件，不仅提高了模拟者的工作效率，也为今后Tough2网格文件的生成提供了新思路和新方法。

4 结论

(1)在ArcGIS平台中利用Python工具箱可在具有地理信息的实体数据上直接生成网格文件，可对任意边界的研究区进行网格剖分，并实现网格的编辑和分析功能。

(2)本次编程是通过生成随机点再生成泰森多边形，对于较为复杂的边界可以利用Python语言将局部网格进行加密，提高模拟的精度。对于其他网格类型，ArcGIS也提供了渔网函数实现矩形网格剖分。

(3)本次编写的Python工具箱，主要针对二维多边形，并以泰森多边形网格的生成为主，可满足众多数字模拟项目的需要。针对三维立体研究区，需借助ArcGIS的三维平台实现二次开发，有待进一步研究。