基于ArcGIS二次開發(fā)的Tough2多邊形網(wǎng)格文件生成方法研究

張 青，胡 成

(中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 4300741)

地下熱流及水流數(shù)值模擬程序Tough2發(fā)布于1991年，經(jīng)過20多年的開發(fā)，已發(fā)展為現(xiàn)今的新流體模塊，應(yīng)用范圍極其廣闊，如CO2地質(zhì)封存、地下水流模擬、環(huán)境影響評價、地?zé)醿Σ毓こ痰?。雖然Tough2可以模擬各類情況，但是推廣上卻較為困難；需要用戶具備一定的編程基礎(chǔ)[1- 3]。在模擬過程的初期，模型構(gòu)建和網(wǎng)格剖分工作量巨大，后人開發(fā)出相應(yīng)的網(wǎng)格生成工具，如MeshMaker、PetraSim、MView等，但是這些網(wǎng)格生成工具并不具有地理信息的概念，也缺乏對網(wǎng)格圖形強大的處理能力。

ArcGIS是具有強大空間分析與處理功能的可視化軟件，可以對具有空間信息的多邊形進行編輯，其平臺提供的二次開發(fā)Python語言是一種跨平臺的開源腳本語言，具有簡單易學(xué)、開發(fā)效率高、伸縮性和可移植性高、可嵌入和穩(wěn)定成熟等優(yōu)勢[4- 5]。自ArcGIS 9.0引入Python作為腳本語言發(fā)展至今，Python已成為一種運用到ArcGIS中實現(xiàn)其數(shù)據(jù)管理與分析等多種功能自動化語言，Python可以創(chuàng)建在獨立腳本環(huán)境、工具箱或模型中運行的腳本。同時，Python可以輕松調(diào)用ArcGIS中各種工具，或訪問其他人員開發(fā)的Python模塊[6- 7]，從而實現(xiàn)程序化處理復(fù)雜數(shù)據(jù)，以減少工作量和提高工作效率。本文根據(jù)數(shù)值模擬項目實際需要，在ArcGIS平臺通過Python實現(xiàn)Tough2多邊形網(wǎng)格的生成，為Tough2的進一步模擬提供便捷的基礎(chǔ)。

1 方法與原理

1.1 Tough2網(wǎng)格文件規(guī)則

Tough2模擬器需要的網(wǎng)格文件包括兩部分，ELEME單元文件和CONNE連接文件，前者反映每個網(wǎng)格單元的基本信息，如單元標(biāo)識、XYZ坐標(biāo)、參數(shù)分區(qū)標(biāo)識等；后者反映每個單元網(wǎng)格之間的關(guān)系，如交界面的面積、網(wǎng)格中心到交界面的距離、重力方向與網(wǎng)格中心連線的余弦值、熱交換因子等，具體參數(shù)見表1和表2。其實，大部分參數(shù)都是網(wǎng)格固有的屬性，但有些參數(shù)是需要根據(jù)實際情況進行賦值的，如參數(shù)分區(qū)標(biāo)識和熱交換因子。實際上Tough2支持任意形狀的網(wǎng)格剖分，無論是三角性、矩形還是任意多邊形都是可以的。

1.2 泰森多邊形生成原理

在一定范圍內(nèi)生成離散點，利用離散點構(gòu)建三角形網(wǎng)，并作出每個三角形各邊的垂直平分線，每個離散點被各邊垂線形成的一個多邊形包圍，這個多邊形即為泰森多邊形[8- 9]。泰森多邊形網(wǎng)格非常符合Tough2模擬器的特點，其不規(guī)則三角網(wǎng)的建立是根據(jù)Delaunay準(zhǔn)則生成，在網(wǎng)格剖分上更具零活性，ArcGIS平臺中自帶泰森多邊形生成工具，利用Python語言可輕松調(diào)用該工具，并對生成的多邊形進行信息提取和圖形分析，如圖1所示。

表1 ELEME單元文件參數(shù)說明

表2 CONNE連接文件參數(shù)說明

圖1 ArcGIS隨機點生成泰森多邊形

1.3 網(wǎng)格文件信息獲取

生成泰森多邊形后可利用ArcGIS提供的內(nèi)置函數(shù)提取網(wǎng)格文件所需要的坐標(biāo)、體積、面積等參數(shù)，并可以利用Python基本函數(shù)計算距離、余弦值等參數(shù)，再根據(jù)實際情況賦予屬性參數(shù)，最后可將提取和計算出的參數(shù)寫入網(wǎng)格文件中。

2 技術(shù)流程

ELEME單元文件和CONNE連接文件正是Tough2所需的MESH網(wǎng)格文件，先有ELEME單元文件，再有CONNE連接文件；即先有網(wǎng)格單元，再有網(wǎng)格關(guān)系。本次按先獲取ELEME文件再獲取CONNE文件進行研究，技術(shù)流程如圖2所示。

圖2 技術(shù)流程圖

2.1 ELEME文件參數(shù)獲取

(1)隨機點的生成

在網(wǎng)格邊界內(nèi)生成一些隨機點，通過隨機點劃分出符合Delaunay準(zhǔn)則的不規(guī)則三角網(wǎng)(TIN)，三角形各邊的垂直平分線即可形成泰森多邊形的邊。在ArcGIS中提供了2種方法用于生成隨機點，一種是使用創(chuàng)建隨機點的工具，可以按照工具給出的提示生成指定范圍內(nèi)的隨機點；另一種是調(diào)用ArcGIS中CreateRandomPoints_management()函數(shù)。在該函數(shù)中只需要給定邊界文件、隨機點個數(shù)、最小點間距等參數(shù)便可快速生成隨機點，其函數(shù)的部分參數(shù)說明見表3。

表3 隨機點函數(shù)參數(shù)說明

生成的隨機點文件中，系統(tǒng)已經(jīng)默認賦予點文件一些屬性，如ID編號、X坐標(biāo)、Y坐標(biāo)等，但這些系統(tǒng)分配好的參數(shù)并不能滿足Tough2網(wǎng)格文件生成所需要的信息，還需在隨機點文件中添加材質(zhì)名稱這一屬性。

(2)創(chuàng)建泰森多邊形

創(chuàng)建泰森多邊形的前提是在已知點位置的條件下，根據(jù)Delaunay準(zhǔn)則生成的。與創(chuàng)建隨機點類似，ArcGIS也提供了2種方法用于創(chuàng)建泰森多邊形：一種是利用CreateThiessenPolygons工具箱；另一種是用CreateThiessenPolygons_analysis()函數(shù)，同樣該函數(shù)只需要輸入隨機點文件和輸出位置便可以生成泰森多邊形，具體參數(shù)見表4。

表4 泰森多邊形函數(shù)參數(shù)說明

與隨機點函數(shù)不同的是，生成的泰森多邊形不需要指定輸出范圍參數(shù)，所以，生成的泰森多邊形是一個矩形，該矩形的范圍是輸入點輸入要素的范圍另加10%，在生成后還需用Clip工具剪裁至模擬區(qū)邊界形狀。

(3)泰森多邊形網(wǎng)格信息提取

生成符合要求的泰森多邊形后需要提取每個單元格的體積，此次研究二維網(wǎng)格的生成，實際上單元體積在數(shù)值上等于單元面積乘以模擬層厚度。X、Y、Z的坐標(biāo)值則需要根據(jù)重力方向進行選取，默認重力方向在Z坐標(biāo)軸上，豎直向下，X坐標(biāo)和Y坐標(biāo)等于投影坐標(biāo)系減去投影帶號所在范圍的基數(shù)值，Z坐標(biāo)等于模擬層厚度的一半。如果重力方向為Y方向，則Y坐標(biāo)與Z坐標(biāo)互換，Y坐標(biāo)在數(shù)值上等于模擬層厚度的一半。對于平行多層網(wǎng)格，在空間上改變的只需是Z值。

在ArcGIS中，單元格面積和坐標(biāo)是圖形固有屬性，上述網(wǎng)格信息可以通過數(shù)據(jù)訪問模塊(arcpy.da)中的Search Cursor類進行只讀權(quán)限訪問。

至此，ELEME文件的主要參數(shù)ID索引碼、X坐標(biāo)、Y坐標(biāo)、單元體積、參數(shù)獲取完成。

2.2 CONNE文件參數(shù)獲取

(1)創(chuàng)建鄰近關(guān)系表

ArcGIS中的面領(lǐng)域工具用于創(chuàng)建領(lǐng)近關(guān)系表，面鄰域工具按照等級路徑確定要在輸出表中記錄的鄰域類型和統(tǒng)計數(shù)據(jù)。鄰域關(guān)系、重合邊和結(jié)點鄰域按照從高到低的等級順序進行重疊。一旦找到更高順序的鄰域，該工具計算和存儲關(guān)系信息并跳過較低順序關(guān)系的分析。由于在創(chuàng)建的泰森多邊形網(wǎng)格(TOUGH2網(wǎng)格)中，領(lǐng)域類型不存在重疊面，需要的是重合邊以獲取相交面積這一參數(shù)。

在圖3中，左邊輸入要素是一個3×3的網(wǎng)格單元，單元編號從101～109，右邊是通過面領(lǐng)域工具分析后，輸出的面領(lǐng)域分析表。

圖3 面領(lǐng)域關(guān)系表

在圖3輸出的面領(lǐng)域關(guān)系表中，揭示了相鄰網(wǎng)格單元間的關(guān)系，與101網(wǎng)格相鄰的有102、104、105網(wǎng)格單元，102和104網(wǎng)格單元的領(lǐng)域類型是重合邊，105網(wǎng)格單元的領(lǐng)域類型是結(jié)點領(lǐng)域。對于Tough2網(wǎng)格文件需要獲取的是重合邊，在圖5中是LENGTH屬性字段。輸出領(lǐng)近關(guān)系表后通過數(shù)據(jù)訪問模塊(arcpy.da)中的SearchCursor訪問這些信息。

(2)獲取公共邊到節(jié)點的距離

每兩個相鄰單元格都只有兩個共同的端點，通過數(shù)據(jù)訪問模塊進行訪問，比較每個網(wǎng)格端點的X坐標(biāo)和Y坐標(biāo)是否相等來獲取公共端點的X坐標(biāo)和Y坐標(biāo)，這樣每個單元格中心與兩個端點間便構(gòu)成了一個三角形，可通過海倫公式：

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

和面積公式S=1/2ah求出單元格中心到公共邊的距離d，即網(wǎng)格中心到公共邊所形成等腰三角形上的高h。

實際中泰森多邊形具有相鄰關(guān)系的兩個單元格中心到公共邊的距離相等，通過兩中心的坐標(biāo)求出距離，其距離的一半即公共邊到節(jié)點的距離。

(3)獲取cosine值

在二維網(wǎng)格中如果取重力方向在Y坐標(biāo)軸上，獲取重力方向與連接線夾角的cosine值方法與計算公共邊到節(jié)點的距離方式類型，以相鄰網(wǎng)格單元的中心節(jié)點(隨機點)構(gòu)建直角三角形，cosine值等于兩中心節(jié)點間的高差比上兩中心節(jié)點的距離。在構(gòu)建直角三角形前需要判斷兩種情況的出現(xiàn)，如果兩節(jié)點的X坐標(biāo)相等，則cosine值等于0，如果Y坐標(biāo)相等，則cosine值等于1或-1。重力方向取X方向與重力方向取Y坐標(biāo)軸方向相似，在此不做說明。如果重力方向在Z坐標(biāo)軸上，則所有的cosine值取0。

至此，CONNE文件的主要參數(shù)相鄰單元格的ID索引碼、單元節(jié)點距離、cosine值的獲取完成，其中，滲透性編號與網(wǎng)格單元的空間位置關(guān)系有關(guān)，如果兩中心節(jié)點在X平面上則為1，在Y平面上則為2，在Z平面上則為3，其反映的是XYZ方向的滲透率，1、2、3只是用于識別滲透性方向的代號，可以在網(wǎng)格文件生成時寫入。

2.3 MESH文件參數(shù)整合

在獲取Tough2中的ELEME和CONNE所需的主要參數(shù)后，將這些參數(shù)按照指定的格式寫入到MESH文件中，使用Python中的format函數(shù)格式化輸出，例如體積的格式化輸出:Volunme=format(ELEME[i][2]，′.4E′)。使用Python中的rjust可以將字符進行右對齊，如:f.write(str(D1).rjust(10))。為了將上述繁瑣的步驟簡單化，方便Tough2網(wǎng)格文件的生成，以Python工具箱的形式進行集成，其工具箱界面如圖4和圖5所示。

圖4 MESH Maker-Polygon.pyt工具箱

圖5 MESH FIle-Polygon.pyt工具箱

圖6的Python工具箱用于生成指定范圍內(nèi)的泰森多邊形網(wǎng)格(Tough2網(wǎng)格)，圖7的Python工具箱用于輸出MESH文件。

2.4 參數(shù)分區(qū)

由于在實際模擬項目中，研究區(qū)內(nèi)各個網(wǎng)格單元的參數(shù)是不同的，在進行模型識別時，不可能也沒有必要讓每個單元都賦予不同的參數(shù)值，所以只能進行參數(shù)分區(qū)，在每個參數(shù)分區(qū)內(nèi)的各單元的參數(shù)是相同的[10]。先在ArcGIS平臺將研究區(qū)劃分為多個子區(qū)域，每個子區(qū)域內(nèi)的單元網(wǎng)格參數(shù)一致，在ELEME單元文件中，每一行有一個位置記錄該單元的參數(shù)，用相應(yīng)的代號表示，可以是英文字母也可以是數(shù)字或者組合，對于具有相同參數(shù)的單元該代號相同。在ArcGIS平臺中可利用Python語言將這些落在同一個參數(shù)分區(qū)內(nèi)的單元節(jié)點的參數(shù)字段統(tǒng)一賦予相同的代號，并對同一參數(shù)分區(qū)內(nèi)的網(wǎng)格進行加密處理。

3 方法應(yīng)用

以安徽某一石油洞庫模擬為例，為了預(yù)測在庫區(qū)施工及運營條件下，庫區(qū)及周邊地下水水位下降情況及其對地質(zhì)環(huán)境可能的影響，需要在分析研究區(qū)水文地質(zhì)條件的基礎(chǔ)上，建立水文地質(zhì)概念模型。根據(jù)地質(zhì)時代、成因、巖性及工程性質(zhì)的不同，將研究區(qū)的巖性分為三大類。其中研究區(qū)東南部主要為第四系地層，成分以粘土、砂土、礫石為主；中部主要為印支期花崗巖體，成分以鉀長石、斜長石、石英、黑云母為主；西北部主要為寒武紀(jì)片巖段，成分以云母石英片巖，頂部為千枚巖為主；洞庫建立在中部的花崗巖體區(qū)域。研究區(qū)參數(shù)分區(qū)如圖6所示。

圖6 研究區(qū)參數(shù)分區(qū)圖

在ArcGIS中勾出模擬區(qū)邊界的矢量文件后，在Python腳本中輸入需要生成的網(wǎng)格數(shù)以及網(wǎng)格中心最小間距即可將模擬區(qū)剖分為泰森多邊形網(wǎng)格，如圖7所示。

圖7 研究區(qū)生成的泰森多邊形

根據(jù)落入每個參數(shù)分區(qū)內(nèi)節(jié)點賦予巖性代號，并自動生成Tough2需要的MESH網(wǎng)格文件。本次處理的坐標(biāo)系為“1954北京坐標(biāo)系”，網(wǎng)格中心最小間距為10m，網(wǎng)格數(shù)量為100個，網(wǎng)格層厚為20m。網(wǎng)格文件結(jié)果如圖8—9所示。

圖8 生成的ELEME網(wǎng)格文件

圖9 生成的CONNE網(wǎng)格文件

在圖8的ELEME單元文件中，第一列代表的數(shù)值含義與表1中第一列的含義一致，即圖8中的數(shù)字1代表單元標(biāo)識為1的網(wǎng)格單元；在圖9的CONNE網(wǎng)格文件中，第一列第一行表示相鄰網(wǎng)格單元1和2的滲透性、距離等信息，分別與表2中的說明一一對應(yīng)。通過在ArcGIS平臺上生成Tough2所需的MESH網(wǎng)格文件，不僅提高了模擬者的工作效率，也為今后Tough2網(wǎng)格文件的生成提供了新思路和新方法。

4 結(jié)論

(1)在ArcGIS平臺中利用Python工具箱可在具有地理信息的實體數(shù)據(jù)上直接生成網(wǎng)格文件，可對任意邊界的研究區(qū)進行網(wǎng)格剖分，并實現(xiàn)網(wǎng)格的編輯和分析功能。

(2)本次編程是通過生成隨機點再生成泰森多邊形，對于較為復(fù)雜的邊界可以利用Python語言將局部網(wǎng)格進行加密，提高模擬的精度。對于其他網(wǎng)格類型，ArcGIS也提供了漁網(wǎng)函數(shù)實現(xiàn)矩形網(wǎng)格剖分。

(3)本次編寫的Python工具箱，主要針對二維多邊形，并以泰森多邊形網(wǎng)格的生成為主，可滿足眾多數(shù)字模擬項目的需要。針對三維立體研究區(qū)，需借助ArcGIS的三維平臺實現(xiàn)二次開發(fā)，有待進一步研究。