地下水數值模擬與可視化建模系統框架研究

閆靜 李立國

摘要：針對現有孔隙地下水有限元數值模擬在三維可視化與空間分析等方面存在的不足，以有限元數值計算方法和3D GIS平臺為基礎，結合GIS空間分析算法和計算機圖形學理論，提出了孔隙地下水有限元分析過程的概念模型的構建、含水層介質的空間離散、水文地質參數提取與賦值等關健步驟，給出了孔隙地下水有限元數值模擬過程在3D GIS下的實現方法和技術框架。該技術框架旨在簡化有限元分析流程、優化模型計算效率、實現地下水有限元數值模擬過程及計算結果的三維可視化。結果表明：基于3DGIS構建的三維水文地質模型，可充分反映地下水系統在3D空間上的分布特征，形成合理的水文地質概念模型，所提出的自適應格網動態生成機制，能夠有效改進有限元數值模擬精度和效率。

關鍵詞：有限元數值模擬；3D GIS；水文地質建模；模型計算參數；可視化

中圖分類號：TV211.1+2 文獻標志碼：A doi：10.3969/j.issn.1000-1379.2018.03.018

地下水有限元數值模擬是對地下水進行定量評估與預測的重要手段，目前國際上已有很多成熟的地下水數值模擬軟件，如MODFLOW、FEFLOW、GMS等，這些軟件功能齊全，能夠有效支持地下水的數值模擬。近年來，隨著GIS理論逐漸成熟與普及，一些學者開始研究利用GIS技術實現地下水數值模擬，Venkatesh Uddameri等建立了多重標準地下水流最優路徑模擬的GIS輔助決策體系；Elkadi A I等基于GIS建立了面向區域特征的地下水流模擬模型；孫繼成等綜合運用FEFLOW軟件和GIS技術，對甘肅省秦王川盆地的水文地質條件進行概化，建立了該地區的地下水數值模擬模型；陳鎖忠等研究了基于GIS的不規則六面體含水層三維空間離散方法，給出了GIS下孔隙地下水流非穩定運動過程數值模型參數的自動提取方法，有效提高了模擬精度。隨著3D GIS技術的發展，數值計算過程及模擬結果的可視化逐漸成為研究熱點，武強等研究了地下水滲流場的動態模擬可視化方法，畢振波等開發了地下水有限元計算可視化系統。

將GIS技術引入數值模擬過程，有利于兼顧數值模型參數的空間特性與精度，提高模型空間分析能力和數據管理效率，使地下水數值模擬更具空間輔助決策特性。筆者基于地下水有限元數值模擬的基本原理，研究了3D GIS下孔隙承壓非穩定地下水有限元數值模擬實現及可視化方法。基于改進前沿推進法（AFT）生成三角形空間離散格網，該方法具有邊界擬合度高、單元尺寸可隨水力梯度自適應的特點；根據三維地質建模理論和空間分析算法，構建三維水文地質模型及其空間分析方法；采用格網約束局部拓撲重建算法提高參數提取精度，基于局部分塊插值技術和計算機圖形學理論實現模型參數的自動賦值；通過三維水流場、三維水文地質模型及其集成形式對計算結果進行可視化。

1 研究區地下水概念模型構建

1.1 研究區水文地質條件

研究區屬于鹽城濱海平原水文地質分區，位于黃海之濱，介于32°75′-34°42′N，119°34′-120°41′E。該平原由2000～3000a來海水和河流不斷沖積而成，區域內水網密布，形成濱海水網平原地貌類型，至今仍在不斷向海延伸。研究區地勢低平，從東南向西北緩慢傾斜。大豐境內地勢較高，海拔3.00～5.00m，向北逐漸降低，到射陽河處為1.00～1.50m。研究區承壓地下水主要賦存于松散沉積巖層孔隙中，埋深為50.00～350.00m。

研究區內含水層多為由粉砂、細砂、粗砂、礫石等組成的松散沉積巖類含水層，含水量充沛，水量穩定，不同含水層之間由黏土或淤泥質黏土層形成隔水層。根據沉積物時代、成因、地層結構以及水文地質特征分析，區內包含孔隙潛水含水層、第Ⅰ承壓含水層、第Ⅱ承壓含水層、第Ⅲ承壓含水層、第Ⅳ承壓含水層5個含水層組（見圖1），其中：孔隙潛水或微承壓水主要接受大氣降水、地表水、農業灌溉水入滲補給，排泄為蒸發和開采利用；第Ⅰ承壓水、第Ⅱ承壓水為微咸水、半咸水、咸水；第Ⅲ承壓含水層由下更新統中細砂、中粗砂組成，厚度為20.00～35.00m，富水性較好，為研究區地下水的主開采層；下伏新近系地層中發育有第Ⅳ、第Ⅴ承壓含水層組，接受上部越流補給和西中部山體的側向補給，主要排泄方式為人工開采。

1.2 地下水連續性運動方程

通常可采用如下連續性運動方程描述地下水的非穩定運動：式中：Kz為越流補給系數；dz為越流經過的垂直距離；Th為含水層水平方向的導水系數；W為單位時間、單位面積上含水層垂向補給水量（流出以負值表示）；S為含水層貯水系數；Hz為水頭；H為待計算水頭；D為滲流區，為一類邊界Г1和二類邊界Г2包圍的研究區域；H（x，y，0）為初始時刻水頭的二維空間分布函數；H0（x，y）為任意位置處初始水頭；q為邊界上單位寬度內流入的流量（流出取負值）；n為外法線方向；t為模擬時段；φ（x，y，t）為指定時間t位置（x，y）處的地下水水位；T為導水系數。

1.3 系統技術框架概述

由上述地下水連續性運動方程數值計算過程可知，3D GIS技術與地下水數值模擬進行結合的切入點主要體現在如下幾個方面。

（1）水文地質概念模型建立。孔隙地下水系統空間分布連續，不同含水層之間的水力聯系受含水層空間結構特征影響明顯，為反映地下水系統在三維空間的分布特征，引入三維地質建模理論輔助構建水文地質概念模型。

（2）空間離散格網生成。離散格網生成是進行地下水有限元數值模擬前處理的重要環節，格網單元的形態、尺寸對有限元模擬精度和計算效率有較大影響，為兼顧計算精度和效率，一般要求三角單元尺寸與不同水力梯度相適應，且相鄰單元之間尺寸過渡盡量平緩，避免出現狹長三角形，同時需顧及含水層空間分布特征以及越流補給等約束條件。筆者采用改進前沿推進法生成模擬區域自適應三角形空間離散格網，綜合考慮水力梯度、邊界條件、含水層結構等約束條件，獲取局部區域最優單元尺寸，動態生成新的自適應空間離散格網。

（3）模型參數提取與初始條件賦值。地下水數值模擬模型參數具有明顯空間分布特性，依據其空間分布特征，抽象為點、線、面矢量數據，對計算單元和水文地質參數矢量數據進行疊置分析，實現計算參數的自動提取。根據模擬初始時刻目標含水層初始水頭空間分布、邊界條件、地下水的開采補給排泄、各水文地質參數空間分布情況，結合空間插值算法將模型計算參數自動賦值到計算單元或格網節點上。

（4）模擬結果可視化。采用2.5D水位DEM和3D地下水流場等可視化表達形式，將地下水流場與三維空間數據模型集成即可得到三維水文地質空間數據模型。

三維空間數據模型由一系列形態各異的體元組成，這些離散的體元一方面能夠擬合水文地質對象的空間結構特征，另一方面是地下水非穩定運動有限元模型與空間數據模型實現耦合的橋梁。這種耦合機制包含如下兩個方面：在有限元數值模擬過程中將計算格網的節點置于空間數據模型所使用體元的節點處，每個節點的參數值可依據相關算法從空間數據模型中自動提取，實現空間數據模型參數向有限元模型的傳輸；隨著模擬時段的推進，數值模擬計算結果不斷發生改變，將計算結果實時回饋到空間數據模型的體元節點上，并驅動其空間結構做出相應調整。通過以上機制，即可實現空間數據模型和地下水數值模型之間參數的雙向傳輸及互相作用。

2 關鍵技術研究

2.1 基于廣義三棱柱（GTP）的三維水文地質建模與空間分析

基于GTP的三維水文地質建模的關鍵步驟：①地質鉆孔建模，整理水文地質鉆孔數據，結合數據特征提取隔水層與含水層信息，構建控制性水文地質鉆孔模型；②含水層界面建模，確定模擬區域邊界范圍，基于空間插值算法構建各含水層頂、底板三角網（TIN）模型；③GTP體元構建，將相鄰含水層界面上的格網節點依據空間關系組織成GTP體元，使得單個GTP體元不與其他GTP體元交叉、重疊，且不會跨越不同類型含水層；④構建含水層模型，GTP體元包含屬性信息，屬性相同的GTP體元組成三維地質體模型。研究區地下水系統及第Ⅲ承壓含水層三維模型見圖2。

為輔助概念模型的空間認知，可進一步開發空間分析模塊，即基于GTP的剖切算法對三維模型進行空間剖切、剖面圖提取、空間開挖等，實現對地下水系統的多角度可視化表達；基于向量場可視化方法，繪制不同時刻的三維地下水流場模型，以反映地下水的時空動態變化特征。

2.2 基于改進前沿推進法（AFT）的自適應空間離散格網構建

AFT應用較為廣泛，所得格網對復雜邊界擬合度高、自適應性好，單元之間尺寸過渡平滑，能夠滿足有限元數值模擬要求。在AFT的格網生成算法中，根據背景格網動態獲取單元尺寸，是實現格網自適應性的關鍵。將水位DEM作為背景格網，建立水力梯度和單元控制尺寸之間的映射關系，預估剖分最大單元尺寸hmax及最小單元尺寸hmin，計算區域內部最大水力梯度Vmax和最小水力梯度Vmin。水力梯度Vi和對應單元尺寸hi之間的關系可用下式近似表示：

水力梯度的計算方法可參閱文獻，在此不再贅述。將式（4）對應關系代入AFT算法，即可得到格網單元尺寸隨水力梯度自適應的空間離散格網。對于地下水數值模擬問題，考慮到區域離散受區域內點、線以及面狀地理要素的約束，且前沿邊推進過程中不可避免地出現“交匯”效應，導致所得離散格網在局部區域可能出現單元尺寸過渡不平滑、三角形形態差等問題，因此引入Laplace光順算子對AFT離散格網進行優化。選取2014年10月第111承壓含水層作為地下水水位 DEM背景格網構建自適應有限元離散格網，預設最大單元尺寸hmax為4km，最小單元尺寸hmin為2km，格網生成結果見圖3。

2.3 基于格網拓撲重建的面狀參數提取

由于水文地質含水層在空間上是連續的，大多地下水數值模擬計算參數在空間上呈面狀分布，因此參照已建概念模型和抽水試驗數據，將模擬含水層劃分為若干個水文地質參數分區，使得每個分區內部含水層近似為均質。參數分區在空間上為不規則多邊形，通過將分區多邊形與離散格網進行疊置分析實現計算參數的自動提取，但存在一個單元同時分布于多個水文地質參數分區的情況，見圖4，單元ABC同時跨越了參數分區P1、P2、P3，此時無法確定單元屬于哪一個分區，給單元參數提取帶來一定困擾。

為保證單元內部水文地質參數的一元性，提出參數分區約束的三角網局部拓撲重建算法，將原空間離散格網轉換為分區邊界線約束的三角格網，使每個單元都完全位于單個參數分區內部。分區邊界線約束格網拓撲重建過程：①對于位于分區邊界線上的一個線段，搜索該線段起點和終點所在的三角形索引；②確定該線段影響到的三角形組成的影響域；③刪除影響域內部的原有三角形，以該線段為起始邊，調用Delaunay算法，重剖分影響域的三角網；④獲取分區邊界線上的下一條線段，重復步驟①一③，直至遍歷所有線段，見圖5。

在拓撲重建得到的約束格網中，每個三角單元僅分布于單一的參數分區內，格網單元內部水文地質條件及參數等有較好的代表性和一元性。此時，采用簡單的射線法即可快速獲取三角形疊合的參數分區，從而完成單元水文地質參數的自動提取。格網拓撲重建前后滲透系數提取結果對比見圖6，可以看出，拓撲重建后格網的計算參數提取結果與參數分區在范圍上更匹配。

2.4 點狀分布的計算參數可視化賦值

對于線狀分布的計算參數，如邊界流量、河流水位、河床滲透系數等，可將其對應的線狀要素通過線約束格網拓撲重建的方式歸人離散格網，過程較簡單，因此重點討論點狀分布參數的可視化賦值方法。

地下水數值模擬模型點狀分布參數主要包括地下水開采量和地下水水位，地下水開采量由模擬區內的開采井決定，地下水水位數據來源于模擬區域內的動態監測井。在地下水流有限元數值計算中，開采量和地下水水位在空間上雖然均以點狀分布，但二者之間的空間分布特征并不一致：地下水水位由含水層的空間分布決定，孔隙地下水含水層具有空間連續性，地下水水位數據在空間上呈連續分布；地下水開采則由人類活動引起，不同開采井之間相互作用有限，開采井在空間上可視作以離散點的形式存在。故對于以上兩種點狀參數，應采用不同的可視化賦值方法。

（1）初始水位賦值。水位監測井可提供水位的實際數據作為預報結果精度評價的依據，在進行離散格網生成之后，需把監測井作為新的格網節點插入離散格網中。如圖7所示，點P表示監測井，獲取包含尸點的單元BCE，搜尋與BCE相鄰的三角形單元，判斷P點是否位于相鄰三角形的外接圓內。若點P位于相鄰三角形外接圓內，則刪除該鄰邊，分別連接P點和相鄰三角形的三個頂點構建新的三角形；若點P不位于相鄰三角形外接圓內，直接連接P與鄰邊的兩個端點構建新三角形。之后根據設定的初始時刻，查詢數據庫中每個監測井在該時刻下的實測水位。調用空間插值算法，形成初始水位場，分配到離散格網節點上，完成初始水位賦值。

（2）開采量賦值。數值模擬中開采量指某單元單位面積上的地下水開采量，開采量的賦值方法：在離散格網中查詢包含開采井的三角單元，計算該單元上單位面積、單位時間的地下水開采量，并將該值賦予該三角單元，其他不包含開采井的單元的開采量賦值為0。

在3D GIS環境下，鹽城市第Ⅲ承壓含水層2014年10月地下水水位及開采量可視化賦值效果見圖8，圖中柱體表示開采井模型，計算格網中紅色區域表示地下水開采量為0，其他顏色則對應不同的開采量，其顏色特征為單調離散著色；水位賦值結果以水位DEM的形式可視化，其顏色特點為連續光滑過渡。

上述可視化過程基于3D GIS平臺，實現了對水文地質參數數據的組織、管理、提取與賦值。不同的賦值模塊以獨立子程序包的形式分別設計開發，封裝成獨立插件集成到數值模擬平臺，實現各模塊之間的松耦合和模塊化，為數值模擬的試算、擬合等提供科學計算可視化工具。

2.5 水文地質參數及其分布特征

為實現地下水流動態變化特征的表達，分析數值模擬以及系統動態建模時有限元數值計算的參數，見表1。

2.6 地下水流速的三維可視化

地下水流模擬輸出結果為離散格網上每個計算節點的水位，一般采用三維流場的形式對地下水流模擬結果進行可視化。根據地下水滲流特點，三維流場的數據源是地下水滲流速度，其方向代表地下水流動方向。對于二維地下水流有限元數值模擬，每個格網節點上的流速可通過X、Y方向上的速度分量：vi、vj合成得到。根據達西定律，計算節點上某方向的滲流速度與節點位置處水力梯度成正比。式中：v為滲流速度；H為水位場；l為某個滲流方向；kl為沿l方向的滲透系數。

求滲流速度的關鍵是獲取表示水位空間分布的函數H的表達式。對于某個節點i，將其周圍一定范圍的水位曲面視作一個二次曲面，該曲面表達式為

Hi（x，y）=b0+b1x+b2y+b3x2+b4xy+b5y2（6）式中：x、y、z為空間位置坐標在三個坐標軸上的投影；b0、b1、…、b5為系數。

該范圍內至少有6個離散格網點，列出這6個計算節點的誤差方程矩陣：式中：x1、x2、…、x6，y1、y2、…、y6為估測值；Z1、Z2、…、Z6為誤差值。

根據平差理論，求得系數矩陣B的解為

B=（MTPM）-1MTPZ（9）式中：P為權重矩陣。

由此可得H（x，y）表達式，分別對其求X，Y方向偏導，結合節點對應的滲透系數可得到節點i處的地下水滲流速度。所取節點越多，曲面越平滑，擬合效果越好。選取“刺狀體”表示水流向量，以顏色表示流速大小，將向量長度設為定值，得到水流體單元在X和Y方向上的分量，將其導入3D GIS平臺實現可視化。

3 應用實例

3.1 模擬系統開發

集成以上核心算法，開發孔隙地下水系統三維數值模擬與可視化平臺。系統分為數據存儲層、數據訪問層和應用系統層：基于Oracle 11g和ArcGIS SDE引擎實現空間數據和屬性數據的一體化建庫與管理；利用C++語言搭建系統框架，實現數據庫的查詢與更新；將OpenGL作為三維場景的渲染引擎，實現數值模擬的可視化管理。

3.2 模擬系統輸出結果

以江蘇省鹽城市第uI承壓含水層中地下水數值模擬為例進行驗證。根據收集到的水位數據以“月”為單位的特點，設模擬起始時間為2014年10月15日，模擬時間步長30d，利用該系統推算4個時段（120d）內的水位變化情況，經過數次參數擬合，代入系統中進行計算，輸出計算結果。研究區15個地下水水位監測井節點的實測與模擬計算結果見表2，可以看出，計算最大誤差為1.96m，最小誤差為0.01m。

3.3 模擬結果精度分析

由表2輸出結果可以看出，在模擬的初始時段（10月），模擬誤差絕對值的平均值為0.013m，第3模擬時段的誤差約為1.213m，隨著時段的推移，模擬誤差呈逐漸增大趨勢。計算水位最大誤差為1.96m，出現在第3個時段。

3.4 預測結果三維可視化

讀取離散格網節點不同時段的水位數據，計算地下水滲流速度，對模擬區域內地下水流場、含水層靜態模型和水位DEM進行集成顯示。不同模擬時段，第Ⅲ承壓含水層地下水水位DEM以及流場在空間的分布情況見圖9、圖10，可以看出，該模型實現了地下水數值模擬計算結果的可視化。

4 結語

結合有限元數值計算原理和承壓地下水含水層水文地質特征，研究了3D GIS環境下有限元數值模擬的關鍵技術與實現方法，并提出了基于3D GIS的地下水流有限元數值模擬技術框架。對現有地下水數值模擬模型的改進主要體現在以下三個方面。

（1）基于三維水文地質建模與分析技術輔助構建水文地質概念模型，可以反映三維孔隙地下水系統在三維空間上的連續分布特征，從而為概念模型構建及其空間分析提供可視化工具。

（2）地下水流有限元計算所用的空間離散格網單元尺寸一般要求能夠與水力梯度相適應，通過建立映射函數，將水力梯度映射到單元尺寸上，并結合AFT算法，提出了隨水力梯度自適應空間離散格網的自動生成機制。應用實例表明該格網能夠自適應反映研究區內水力梯度的空間分布，并有效支持有限元計算。

（3）針對模擬模型計算參數空間分布特征，將其劃分為“點、線、面”三種矢量數據類型，并提出了每種數據類型參數的可視化提取與賦值方法。實現了具有較強空間特性的參數數據的高效組織、管理存儲與可視化，避免了繁瑣參數數據文件的組織與管理工作。

此外，完成了動態自適應的地下水自動模擬系統框架。針對松散沉積含水層空間分布特點，基于GTP空間數據模型實現地下水系統自適應三維動態建模。引入AFT算法對動態建模過程中可能遇到的前沿邊的幾何形態類型進行分類和歸納，提出改進AFT的自適應格網構建方法。在此基礎上梳理了地下水滲流運動的有限元數值求解步驟，研究了利用GIS實現各種類型參數數據管理與組織的方法，解決了地下水滲流數值模擬過程模型參數數據管理不便、效率較低的問題，為三維水文地質動態建模及有限元數據模擬提供了高效、靈活、直觀的數據管理模式與平臺。但仍存在如下不足之處：將離散格網單元尺寸和水力梯度之間關系近似為簡單的線性反比例關系得到二者之間映射關系模型，顧及影響因素較少，有待細化；采用“刺狀體”繪制點箭頭的方式對地下水流空間分布特征進行建模，當采樣點較為密集時，容易引發視覺上的混亂，而采樣不足又可能導致關鍵特征的丟失，難以全面、連續地反映向量場。采用流線法、數值積分法、流函數法等，可對上述問題進行修正，基于不同方法實現地下水流場的三維可視化是后續研究的一個重要方向。