閆靜 李立國

摘要：針對(duì)現(xiàn)有孔隙地下水有限元數(shù)值模擬在三維可視化與空間分析等方面存在的不足，以有限元數(shù)值計(jì)算方法和3D GIS平臺(tái)為基礎(chǔ)，結(jié)合GIS空間分析算法和計(jì)算機(jī)圖形學(xué)理論，提出了孔隙地下水有限元分析過程的概念模型的構(gòu)建、含水層介質(zhì)的空間離散、水文地質(zhì)參數(shù)提取與賦值等關(guān)健步驟，給出了孔隙地下水有限元數(shù)值模擬過程在3D GIS下的實(shí)現(xiàn)方法和技術(shù)框架。該技術(shù)框架旨在簡(jiǎn)化有限元分析流程、優(yōu)化模型計(jì)算效率、實(shí)現(xiàn)地下水有限元數(shù)值模擬過程及計(jì)算結(jié)果的三維可視化。結(jié)果表明：基于3DGIS構(gòu)建的三維水文地質(zhì)模型，可充分反映地下水系統(tǒng)在3D空間上的分布特征，形成合理的水文地質(zhì)概念模型，所提出的自適應(yīng)格網(wǎng)動(dòng)態(tài)生成機(jī)制，能夠有效改進(jìn)有限元數(shù)值模擬精度和效率。

關(guān)鍵詞：有限元數(shù)值模擬；3D GIS；水文地質(zhì)建模；模型計(jì)算參數(shù)；可視化

中圖分類號(hào)：TV211.1+2 文獻(xiàn)標(biāo)志碼：A doi：10.3969/j.issn.1000-1379.2018.03.018

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

將GIS技術(shù)引入數(shù)值模擬過程，有利于兼顧數(shù)值模型參數(shù)的空間特性與精度，提高模型空間分析能力和數(shù)據(jù)管理效率，使地下水?dāng)?shù)值模擬更具空間輔助決策特性。筆者基于地下水有限元數(shù)值模擬的基本原理，研究了3D GIS下孔隙承壓非穩(wěn)定地下水有限元數(shù)值模擬實(shí)現(xiàn)及可視化方法。基于改進(jìn)前沿推進(jìn)法（AFT）生成三角形空間離散格網(wǎng)，該方法具有邊界擬合度高、單元尺寸可隨水力梯度自適應(yīng)的特點(diǎn)；根據(jù)三維地質(zhì)建模理論和空間分析算法，構(gòu)建三維水文地質(zhì)模型及其空間分析方法；采用格網(wǎng)約束局部拓?fù)渲亟ㄋ惴ㄌ岣邊?shù)提取精度，基于局部分塊插值技術(shù)和計(jì)算機(jī)圖形學(xué)理論實(shí)現(xiàn)模型參數(shù)的自動(dòng)賦值；通過三維水流場(chǎng)、三維水文地質(zhì)模型及其集成形式對(duì)計(jì)算結(jié)果進(jìn)行可視化。

1 研究區(qū)地下水概念模型構(gòu)建

1.1 研究區(qū)水文地質(zhì)條件

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

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

1.2 地下水連續(xù)性運(yùn)動(dòng)方程

通?？刹捎萌缦逻B續(xù)性運(yùn)動(dòng)方程描述地下水的非穩(wěn)定運(yùn)動(dòng)：式中：Kz為越流補(bǔ)給系數(shù)；dz為越流經(jīng)過的垂直距離；Th為含水層水平方向的導(dǎo)水系數(shù)；W為單位時(shí)間、單位面積上含水層垂向補(bǔ)給水量（流出以負(fù)值表示）；S為含水層貯水系數(shù)；Hz為水頭；H為待計(jì)算水頭；D為滲流區(qū)，為一類邊界Г1和二類邊界Г2包圍的研究區(qū)域；H（x，y，0）為初始時(shí)刻水頭的二維空間分布函數(shù)；H0（x，y）為任意位置處初始水頭；q為邊界上單位寬度內(nèi)流入的流量（流出取負(fù)值）；n為外法線方向；t為模擬時(shí)段；φ（x，y，t）為指定時(shí)間t位置（x，y）處的地下水水位；T為導(dǎo)水系數(shù)。

1.3 系統(tǒng)技術(shù)框架概述

由上述地下水連續(xù)性運(yùn)動(dòng)方程數(shù)值計(jì)算過程可知，3D GIS技術(shù)與地下水?dāng)?shù)值模擬進(jìn)行結(jié)合的切入點(diǎn)主要體現(xiàn)在如下幾個(gè)方面。

（1）水文地質(zhì)概念模型建立。孔隙地下水系統(tǒng)空間分布連續(xù)，不同含水層之間的水力聯(lián)系受含水層空間結(jié)構(gòu)特征影響明顯，為反映地下水系統(tǒng)在三維空間的分布特征，引入三維地質(zhì)建模理論輔助構(gòu)建水文地質(zhì)概念模型。

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

（3）模型參數(shù)提取與初始條件賦值。地下水?dāng)?shù)值模擬模型參數(shù)具有明顯空間分布特性，依據(jù)其空間分布特征，抽象為點(diǎn)、線、面矢量數(shù)據(jù)，對(duì)計(jì)算單元和水文地質(zhì)參數(shù)矢量數(shù)據(jù)進(jìn)行疊置分析，實(shí)現(xiàn)計(jì)算參數(shù)的自動(dòng)提取。根據(jù)模擬初始時(shí)刻目標(biāo)含水層初始水頭空間分布、邊界條件、地下水的開采補(bǔ)給排泄、各水文地質(zhì)參數(shù)空間分布情況，結(jié)合空間插值算法將模型計(jì)算參數(shù)自動(dòng)賦值到計(jì)算單元或格網(wǎng)節(jié)點(diǎn)上。

（4）模擬結(jié)果可視化。采用2.5D水位DEM和3D地下水流場(chǎng)等可視化表達(dá)形式，將地下水流場(chǎng)與三維空間數(shù)據(jù)模型集成即可得到三維水文地質(zhì)空間數(shù)據(jù)模型。

三維空間數(shù)據(jù)模型由一系列形態(tài)各異的體元組成，這些離散的體元一方面能夠擬合水文地質(zhì)對(duì)象的空間結(jié)構(gòu)特征，另一方面是地下水非穩(wěn)定運(yùn)動(dòng)有限元模型與空間數(shù)據(jù)模型實(shí)現(xiàn)耦合的橋梁。這種耦合機(jī)制包含如下兩個(gè)方面：在有限元數(shù)值模擬過程中將計(jì)算格網(wǎng)的節(jié)點(diǎn)置于空間數(shù)據(jù)模型所使用體元的節(jié)點(diǎn)處，每個(gè)節(jié)點(diǎn)的參數(shù)值可依據(jù)相關(guān)算法從空間數(shù)據(jù)模型中自動(dòng)提取，實(shí)現(xiàn)空間數(shù)據(jù)模型參數(shù)向有限元模型的傳輸；隨著模擬時(shí)段的推進(jìn)，數(shù)值模擬計(jì)算結(jié)果不斷發(fā)生改變，將計(jì)算結(jié)果實(shí)時(shí)回饋到空間數(shù)據(jù)模型的體元節(jié)點(diǎn)上，并驅(qū)動(dòng)其空間結(jié)構(gòu)做出相應(yīng)調(diào)整。通過以上機(jī)制，即可實(shí)現(xiàn)空間數(shù)據(jù)模型和地下水?dāng)?shù)值模型之間參數(shù)的雙向傳輸及互相作用。

2 關(guān)鍵技術(shù)研究

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

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

為輔助概念模型的空間認(rèn)知，可進(jìn)一步開發(fā)空間分析模塊，即基于GTP的剖切算法對(duì)三維模型進(jìn)行空間剖切、剖面圖提取、空間開挖等，實(shí)現(xiàn)對(duì)地下水系統(tǒng)的多角度可視化表達(dá)；基于向量場(chǎng)可視化方法，繪制不同時(shí)刻的三維地下水流場(chǎng)模型，以反映地下水的時(shí)空動(dòng)態(tài)變化特征。

2.2 基于改進(jìn)前沿推進(jìn)法（AFT）的自適應(yīng)空間離散格網(wǎng)構(gòu)建

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

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

2.3 基于格網(wǎng)拓?fù)渲亟ǖ拿鏍顓?shù)提取

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

為保證單元內(nèi)部水文地質(zhì)參數(shù)的一元性，提出參數(shù)分區(qū)約束的三角網(wǎng)局部拓?fù)渲亟ㄋ惴?，將原空間離散格網(wǎng)轉(zhuǎn)換為分區(qū)邊界線約束的三角格網(wǎng)，使每個(gè)單元都完全位于單個(gè)參數(shù)分區(qū)內(nèi)部。分區(qū)邊界線約束格網(wǎng)拓?fù)渲亟ㄟ^程：①對(duì)于位于分區(qū)邊界線上的一個(gè)線段，搜索該線段起點(diǎn)和終點(diǎn)所在的三角形索引；②確定該線段影響到的三角形組成的影響域；③刪除影響域內(nèi)部的原有三角形，以該線段為起始邊，調(diào)用Delaunay算法，重剖分影響域的三角網(wǎng)；④獲取分區(qū)邊界線上的下一條線段，重復(fù)步驟①一③，直至遍歷所有線段，見圖5。

在拓?fù)渲亟ǖ玫降募s束格網(wǎng)中，每個(gè)三角單元僅分布于單一的參數(shù)分區(qū)內(nèi)，格網(wǎng)單元內(nèi)部水文地質(zhì)條件及參數(shù)等有較好的代表性和一元性。此時(shí)，采用簡(jiǎn)單的射線法即可快速獲取三角形疊合的參數(shù)分區(qū)，從而完成單元水文地質(zhì)參數(shù)的自動(dòng)提取。格網(wǎng)拓?fù)渲亟ㄇ昂鬂B透系數(shù)提取結(jié)果對(duì)比見圖6，可以看出，拓?fù)渲亟ê蟾窬W(wǎng)的計(jì)算參數(shù)提取結(jié)果與參數(shù)分區(qū)在范圍上更匹配。

2.4 點(diǎn)狀分布的計(jì)算參數(shù)可視化賦值

對(duì)于線狀分布的計(jì)算參數(shù)，如邊界流量、河流水位、河床滲透系數(shù)等，可將其對(duì)應(yīng)的線狀要素通過線約束格網(wǎng)拓?fù)渲亟ǖ姆绞綒w人離散格網(wǎng)，過程較簡(jiǎn)單，因此重點(diǎn)討論點(diǎn)狀分布參數(shù)的可視化賦值方法。

地下水?dāng)?shù)值模擬模型點(diǎn)狀分布參數(shù)主要包括地下水開采量和地下水水位，地下水開采量由模擬區(qū)內(nèi)的開采井決定，地下水水位數(shù)據(jù)來源于模擬區(qū)域內(nèi)的動(dòng)態(tài)監(jiān)測(cè)井。在地下水流有限元數(shù)值計(jì)算中，開采量和地下水水位在空間上雖然均以點(diǎn)狀分布，但二者之間的空間分布特征并不一致：地下水水位由含水層的空間分布決定，孔隙地下水含水層具有空間連續(xù)性，地下水水位數(shù)據(jù)在空間上呈連續(xù)分布；地下水開采則由人類活動(dòng)引起，不同開采井之間相互作用有限，開采井在空間上可視作以離散點(diǎn)的形式存在。故對(duì)于以上兩種點(diǎn)狀參數(shù)，應(yīng)采用不同的可視化賦值方法。

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

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

在3D GIS環(huán)境下，鹽城市第Ⅲ承壓含水層2014年10月地下水水位及開采量可視化賦值效果見圖8，圖中柱體表示開采井模型，計(jì)算格網(wǎng)中紅色區(qū)域表示地下水開采量為0，其他顏色則對(duì)應(yīng)不同的開采量，其顏色特征為單調(diào)離散著色；水位賦值結(jié)果以水位DEM的形式可視化，其顏色特點(diǎn)為連續(xù)光滑過渡。

上述可視化過程基于3D GIS平臺(tái)，實(shí)現(xiàn)了對(duì)水文地質(zhì)參數(shù)數(shù)據(jù)的組織、管理、提取與賦值。不同的賦值模塊以獨(dú)立子程序包的形式分別設(shè)計(jì)開發(fā)，封裝成獨(dú)立插件集成到數(shù)值模擬平臺(tái)，實(shí)現(xiàn)各模塊之間的松耦合和模塊化，為數(shù)值模擬的試算、擬合等提供科學(xué)計(jì)算可視化工具。

2.5 水文地質(zhì)參數(shù)及其分布特征

為實(shí)現(xiàn)地下水流動(dòng)態(tài)變化特征的表達(dá)，分析數(shù)值模擬以及系統(tǒng)動(dòng)態(tài)建模時(shí)有限元數(shù)值計(jì)算的參數(shù)，見表1。

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

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

求滲流速度的關(guān)鍵是獲取表示水位空間分布的函數(shù)H的表達(dá)式。對(duì)于某個(gè)節(jié)點(diǎn)i，將其周圍一定范圍的水位曲面視作一個(gè)二次曲面，該曲面表達(dá)式為

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

該范圍內(nèi)至少有6個(gè)離散格網(wǎng)點(diǎn)，列出這6個(gè)計(jì)算節(jié)點(diǎn)的誤差方程矩陣：式中：x1、x2、…、x6，y1、y2、…、y6為估測(cè)值；Z1、Z2、…、Z6為誤差值。

根據(jù)平差理論，求得系數(shù)矩陣B的解為

B=（MTPM）-1MTPZ（9）式中：P為權(quán)重矩陣。

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

3 應(yīng)用實(shí)例

3.1 模擬系統(tǒng)開發(fā)

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

3.2 模擬系統(tǒng)輸出結(jié)果

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

3.3 模擬結(jié)果精度分析

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

3.4 預(yù)測(cè)結(jié)果三維可視化

讀取離散格網(wǎng)節(jié)點(diǎn)不同時(shí)段的水位數(shù)據(jù)，計(jì)算地下水滲流速度，對(duì)模擬區(qū)域內(nèi)地下水流場(chǎng)、含水層靜態(tài)模型和水位DEM進(jìn)行集成顯示。不同模擬時(shí)段，第Ⅲ承壓含水層地下水水位DEM以及流場(chǎng)在空間的分布情況見圖9、圖10，可以看出，該模型實(shí)現(xiàn)了地下水?dāng)?shù)值模擬計(jì)算結(jié)果的可視化。

4 結(jié)語

結(jié)合有限元數(shù)值計(jì)算原理和承壓地下水含水層水文地質(zhì)特征，研究了3D GIS環(huán)境下有限元數(shù)值模擬的關(guān)鍵技術(shù)與實(shí)現(xiàn)方法，并提出了基于3D GIS的地下水流有限元數(shù)值模擬技術(shù)框架。對(duì)現(xiàn)有地下水?dāng)?shù)值模擬模型的改進(jìn)主要體現(xiàn)在以下三個(gè)方面。

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

（2）地下水流有限元計(jì)算所用的空間離散格網(wǎng)單元尺寸一般要求能夠與水力梯度相適應(yīng)，通過建立映射函數(shù)，將水力梯度映射到單元尺寸上，并結(jié)合AFT算法，提出了隨水力梯度自適應(yīng)空間離散格網(wǎng)的自動(dòng)生成機(jī)制。應(yīng)用實(shí)例表明該格網(wǎng)能夠自適應(yīng)反映研究區(qū)內(nèi)水力梯度的空間分布，并有效支持有限元計(jì)算。

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

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