地下水滲流場可視化方法與應用研究*

何 亮,郜沐晨,陳鎖忠,齊 慧

(1.虛擬地理環境教育部重點實驗室(南京師范大學)，江蘇 南京，210023; 2.江蘇省地理環境演化國家重點實驗室培育建設點，江蘇 南京，210023; 3.江蘇省地理信息資源開發與利用協同創新中心，江蘇 南京，210023)

1 引言

地下水是水資源的重要組成部分，隨著我國人口急劇增長，經濟與城市的快速發展，水資源的需求量急速增加，一些區域與城市地下水開采強度持續增加，使地下水主采層長期處于超采狀態，導致水位大幅度下降，引發了如地下水資源衰減、地面沉降、地裂縫、地面塌陷、海水入侵與水質惡化等系列環境地質問題，嚴重制約了地區經濟的可持續發展[1]。為保護地下水資源與生態環境，各級政府從行政措施與技術手段2個方面加強了地下水合理開采的管理。

地下水滲流場是滲透水流所占有的空間區域，在滲流場內任一點均具有一定的水頭和滲透速度，水頭和滲透速度是滲流區內點的坐標和時間的函數，滲流場具有明顯的時空分異性。基于地下水滲流場的時空演變情況可以分析得到相應區域地下水的開采強度與水位下降程度，是地下水資源合理開發方案制定的重要依據[2]。由于地下水埋藏在地表以下，相比地表水而言，有著更為復雜的形成與運移機理，技術與管理人員不能直接觀測地下水的運動規律、水位動態變化情況，其賦存條件、水位動態變化只能通過水文地質勘查與地下水動態長期監測數據揭示。地下水滲流場傳統的繪制方法是基于水文地質剖面或平面水文地質圖，根據地下水位等值線繪制而成，僅反映了某個時刻的地下水滲流情況，可視化程度較低，不能反映滲流場連續變化的動態特征。

為了增強二維流場中流線的清晰度和可行度，且能反映滲流場連續的時空動態特征，基于滲流場的繪制原理，采用計算機圖形學、GIS空間分析與可視化技術，自動繪制地下水滲流場及其動態演變場景。本文從地下水滲流的流線繪制、流速計算、流場種子點與終止點選擇、流線追蹤4個方面研究孔隙地下水滲流場可視化算法，實現了孔隙地下水二維可視化滲流場的自動構建，同時，提出了基于地下水可視化滲流場的地下水紅線管理方法。

2 孔隙地下水滲流場可視化方法

孔隙地下水滲流場的可視化對地下水動態在空間上表達具有關鍵作用，通過滲流場及其動態演變可分析出地下水流向、流速與水力梯度的空間分異情況。根據場的物理性質，可將場分為向量場與標量場；根據場的物理狀態和時間關系，可將場分為靜態場與動態場[3]。孔隙地下水滲流場屬于向量場。向量場可視化方法包括：點箭頭法、基于紋理法、基于特征法和向量線法等。

點箭頭法作為一種直觀簡單的可視化方法被廣泛應用于流場可視化研究中，其形狀類似于刺狀體。點箭頭法包含3大要素：箭頭、錐體與向量線段，點箭頭法在2D向量場內能夠呈現出滲流區內的流速大小與滲流方向。但是，若點箭頭數量太多或者分布太密集將會使得結果圖中箭頭雜亂無序；若點箭頭數量過少或者分布密度較稀疏，將難以有效呈現向量場的空間分布形態，且數據的內在連續性也無法得到準確表達。所以，點箭頭法僅適用于小區域的流場可視化。

紋理法是用極其精細且連貫的密集紋理線條、顏色繪制與渲染可視化流場，包括點噪聲法和線積分卷積法。紋理法不但可以生成高密度的可視化效果，同時也可以表征全局流場。和傳統的點箭頭法比較，該方法可以更好地呈現其細節特征，在二維流場與曲面流場中有著廣泛應用。但是，由于紋理法可視化效果比較精細，數據計算量非常大，需要更高性能的處理器。

特征法大多用于對某一具體問題的判別，針對某一具體情況，特征可廣義地定義為任一實體、結構或范圍[4-6]，所以大部分基于特征的可視化方法是針對某個特殊問題。特征法可分為基于圖像處理、基于空間拓撲和基于物理特征3類。

向量線法中，用完全正切于向量場的流線段來表征流場。基于向量線的可視化形式分為流線、脈線和跡線。穩定向量場一般選用流線進行可視化，利用數值積分獲取向量曲線來表征二維流場。不同的積分曲線以及種子點的選取方法不同是引起向量線法之間差異的重要因素。相較于點箭頭法，向量線法因為線的連貫性，在表征場的連貫性方面有更好的可視化效果。本文采用向量線法對地下水流場特征進行可視化。

3 流線的繪制

3.1 流線的定義

地下水流線可反映地下水滲流場在時空域的變化情況，流線為滲流場中某個瞬時的一條線，線上各點在此瞬時的流線均與此線相切。按照流線的定義，可得到在流場中，瞬時t0經過點P0的一條流線。首先，在點P0處繪制速度矢量v0，并順著v0速度矢量方向，距點P0長度為▽a的點P1，繼續繪制該時刻的速度矢量v1，接著順著v1速度矢量方向，距點P1長度為▽a的點P2，繪制該時刻的速度矢量v2，再順著v2速度矢量方向，距點P2長度為▽a的點P3，繪制該時刻速度矢量v3，用相同的方式繼續繪制，最終獲得折線P0P1P2P3P4，設2個點的長度▽a逼近于0，則折線P0P1P2P3P4會轉換成一條光滑的曲線，依此繪制的曲線叫做流線，如圖1所示。

Figure 1 Streamline definition圖1 流線的定義

表達流線的微分方程能通過流線定義得出[7]。在二維空間中，點的速度矢量與流線的關系是相切的關系，故速度矢量v和該流線微元弧的矢量積為0。在流線上提取微小線元ds，此時速度方向為微小線元ds的方向，如果2個矢量平行，則其叉乘值為0，即：

v×ds=0

(1)

由于流線是一條曲線，并且與滲流速度矢量相切，流線簇示意滲流區中任意點處的水流運動方向。在任意流線上取任意2點P(x,y)和P′(x+dx,y+dy)。P點的滲流速度矢量為v，它與它的2個分量vx和vy構成一個三角形PMN。自P′作垂線P′m，并延長至n，如圖2所示。

當P(x,y)與P′(x+dx,y+dy)無限接近時，弧線PP′可用切線Pm來替換。即Pn=dx,mn=dy。又因ΔPMN～ΔPmn，所以：

vxdy-vydx=0

(2)

Figure 2 Streamline diagram圖2 流線的示意圖

由于P(x,y)與P′(x+dx,y+dy)是此流線上隨機的2個點，所以針對流線上任何點，式(2)皆有效，即能把它視為流線方程，可用于描述流線[8]。而式(2)僅表征某點的水流運動狀態，故均適用于各向同性、異性介質，結合達西定律廣義公式，則式(2)可變為：

(3)

其中，Kxx和Kyy分別為X方向和Y方向的滲透系數，H為地下水位置。對于各向同性介質，式中Kxx=Kyy。而對于各向異性介質，選取坐標軸(直角坐標系)的方向分別與滲透系數的主方向一致的情況。

3.2 流速的計算

流速計算準確性會直接影響最終流線的準確性。地下水流速(滲透速度)指的是當水流經過含骨架和孔隙中的斷面時存在的一類模擬流速。法國水力學家Darey于1856年發現，孔隙地下水是從水位高處流向水位低處，最后得出結論：

(4)

q=ui+vj+wk

(5)

水力坡度J在x,y,z方向上的分量Jx,Jy,Jz可表示為：

(6)

在梯度上，曲面中任意點的梯度和經過該點的等值線之間存在這樣的聯系：即函數z=f(x,y)于P(x,y)處的梯度方向，和經過P(x,y)處的等高線f(x,y)=C于此點的法線方向一致，遵循從等值線低值指向等值線高值的規律，另外梯度的模為函數于此法線方向上的導數，曲面上一點的梯度與過這點等值線的關系如圖3所示。

Figure 3 Relationship between the gradient of a point on a surface and the contour of this point圖3 曲面上一點的梯度與過這點等值線的關系

根據水力梯度和梯度之間聯系，先統計曲面等值線梯度值，再得出水力梯度值，方法如下所示:

將三維空間的地下水位等值面(曲面)變化特征投影到二維上是一個平面，將二維平面(x,y)坐標作為自變量，地下水位H作為因變量。函數H=f(x,y)反映或趨近其地下水位的地域分布情況，而統計流速所涉及的水位數據必須至少有3個及以上，本文采取局部趨勢面法來擬合水位數據。為改善趨勢面的擬合度，并避免出現所求點附近的現水位點數量偏少的情況，本文構建了能夠涵蓋所有研究范圍的規則格網，再結合監測井的水位值，用樣條函數插值法統計研究范圍中在格網交點處的水位。

通過樣條函數插值統計得到格網點的水位與原始監測井的水位，并將其一并用于建立水位動態變化趨勢面。若所求點是P(x,y)，則定義半徑R，并以P為圓心，以R為半徑畫圓，圓中所有涵蓋的采樣點用于設立的水位趨勢面。考慮全局趨勢與擬合效果，選用次數為2的趨勢面模型進行水位曲面的統計。此時，趨勢面上P(x,y)點的梯度向量是：

(a1+2a3x+a4y,a2+2a4x+a5y)

(7)

又因地下水是從機械能較高點流向機械能較低點，故水力坡度與梯度向量的方向相反：

J=(-a1-2a3x-a4y,-a2-2a4x-a5y)

(8)

之后再按照式(4)計算P點的速度。

3.3 種子點和終止點的選擇

合理的流線種子點與終止點選取策略是流線空間分布、洞悉流場特征的重要影響因素。通過限制初始種子點與終止種子點的空間位置來控制流線的生成，首先構建并生成種子點集合，然后遍歷種子點集合中所有種子點，并按照其流向，利用數值積分的計算方法來獲取流線。只要種下充分數量的種子點，即可在流場中產生充分數量的流線，從而不存在流場特性的丟失情況。然而，如果種子點與終止點的選擇方案不合理，會出現流線無序和大面積留白的問題，直接影響可視化結果。

一般的種子播種方法具有獨特的優勢和劣勢，包括隨機放置法、等距放置法、流場拓撲誘導放置法和同類散播放置法等方法。考慮到布局策略的利弊，本文采用一種新的種子點放置策略來表征地下水流場，主要思路為地下水從高水位處流向低水位處。通過3個因素來確定種子點的選擇：

(1)地下水流動方向。

地下水流場作為一種地下水對機械能差別的呼應，地下水自具有較高機械能點流向具有較低機械能點，通常順著最陡的水力梯度進行移動[7]。以前在缺少科學繪圖工具繪制流線的時候，大部分的水文地質學家都是手動繪制，主要方法是基于等水位線自地下水水位值較高的等值線朝地下水水位值較低的等值線繪制流線，以此來展示地下水的流動方向。

(2)重點研究區。

地下水的超限采地區、供水水源區、生態脆弱地區、水資源治理和保障地區等作為重點研究區，重點研究區的地下水水位變化會直接影響該區域的資源、環境和生態等，所以當對地下水流場進行科學可視化時，需對此重點研究區的水流狀況進行可視化。

(3)流線繪制合理性。

理論上，如果種下充分的種子，就能產生充分的流線，故所有流場的特性就不會損失。但是，若流線的長度太長、播種的種子數量過多都將會導致流線混淆或隱藏，部分可能遮住流場的特性，所以應綜合考慮流場的特性與科學可視化的效果來確定流線的具體條數。通常，當相交的等水位線和流線2者之間構成近似正方形單元時效果最佳。

綜上因素，本文種子點的選取方法如下所示：

(1)初始種子點的選擇。

基于規則的格網與最大等水位線來篩選所需種子點，并通過規則格網的數量來判斷和描述小網格涵蓋范圍內是否存在流線，其中小網格長寬根據所需流線的密度來選擇。為每一個小網格定義其唯一的標識變量，以計算通過小網格的流線數量，且初始值設置為零。

把規則的格網和最大等水位線在空間上進行疊加分析，并將最大等水位線通過的小網格的中心點當做種子點，且該小網格的計數變量加1。在選擇初始種子點之后，即可以開始流線追蹤，生成并繪制流線，在該過程中，對通過流線的小網格的計數變量遞增1。其中初始種子點選擇的示意圖如圖4所示。

Figure 4 Initial seed point selection diagram圖4 初始種子點選擇示意圖

(2)興趣區種子點的選擇。

為避免初始種子點所產生的流線未能將地質工作者想重點了解的區域詳盡顯示等問題，初始種子點選擇好并生成流線之后，必須校對是否存在流線通過興趣區。若為空白，則需在興趣區中播種種子點。通常，水位動態監測井所處的位置是水文地質工作者重點關注地區，亦或是能夠顯示地下水特性的代表范圍。故將地下水水位監測井數據和規則格網數據疊加分析，若無流線通過地下水水位監測井點附近地區，則把此地下水水位監測井點看成興趣區內的種子點，繪制流線。

(3)加密種子點的選擇。

流線生成之后，需對所有小網格進行遍歷，若存在大塊區域的小網格無流線通過情況，而該范圍內又存在等水位線，此時繼續遵循(1)中策略，將該范圍內的最大等水位線和網格進行疊加分析，并選取種子點作為加密種子點，并繪制流線。

因為流線是有起始點和結束點的一條線，故合理化設定流線結束條件也很關鍵。如果達到以下3個條件中的任何1個，必須終止對流線的跟蹤：當發現流線的積分超出了流場范圍的界限；此刻處于積分點的流速是零；相鄰流線的距離低于設定的閾值。

3.4 流線的跟蹤

我們從上述流線定義發現，很難用具體的顯式函數去描述流線。針對此，本文采用數值積分法，依照確認的種子點，通過流場、流速的計算辦法，開始跟蹤計算，能夠用多個離散點來繪制曲線。一般數值積分法包括歐拉(EULER)方法、改進的歐拉方法和四階Runge-Kutta方法[9 - 11]等。歐拉方法是數值計算中最基本、最簡單的方法之一，然而其計算精度不高，通常不單獨運用于工程計算中。后來有學者對歐拉方法進行了優化，即改進的歐拉方法。首先通過歐拉方法得到預報值，然后使用梯形公式計算得到校正值，發現局部截斷的誤差比歐拉方法的小一階，大大提高了結果的準確性。四階Runge-Kutta方法具有較高的計算準確性，而計算復雜度較高且耗時長。最終，考慮到效率和精度等問題，本文選擇改進的歐拉方法計算：

(9)

流線生成具體技術路線圖如圖5所示。

Figure 5 Streamline generation technology roadmap圖5 流線生成技術路線圖

4 滲流場在地下水位紅線管理中的應用

地下水水位紅線是根據數值模擬計算方法，通過計算開采層的禁采水位來確定的。當含水層的水位低于禁采水位紅線后，持續開采將會誘發系列環境地質問題。

基于等水位面和流線對地下水流場進行科學可視化，一方面，等水位面能夠體現地下水位在空間上的分布形態，另一方面，流線箭頭方向可以表明地下水運動趨勢和地下水之間的補給關系。地下水流場可視化方法在地下水位紅線管理中的具體應用方法如下所示：

(1)水位紅線(禁采埋深、限采埋深)可視化分析方法。相同的水位紅線管理研究范圍中，其含水層的紅線指數保持一致，因此水位紅線指標面也可以看成是一個等水位面[9 - 12]。考慮到地下水位紅線的物理意義與幾何特征，本文運用等水位面生成方法對研究范圍進行了可視化分析，并制作了研究范圍內的限采水位紅線等水位面與禁采水位紅線等水位面。此外，基于GIS的空間數據符號化，在等水位面內，把水位值等于限采水位紅線的范圍用淺灰色進行渲染，該淺灰色區域表示水位值小于限采水位紅線的區域；把水位值等于禁采水位紅線的范圍用深灰色進行渲染，該深灰色區域表示水位低于禁采水位紅線的區域，其余為水位合理區域。基于此方法，即可對研究范圍內的地下水位實行紅線分析。

(2)量化和可視化水位紅線分析結果的方法。基于研究范圍的紅線分析結果，可看出具體哪一個范圍內的水位低于水位紅線以及該區域的地理位置。再基于GIS的空間拓撲方法，拓撲分析水位紅線的等值線圖[12,13]。又由于等水位線存在2種形態，即封閉狀態與非封閉狀態，所以拓撲建立水位等值面過程中需優化：

①非封閉的等水位線：按照拓撲數據結構，一個面是由封閉的線段組成，另外，因為邊界的存在會引發部分等值線沒有封閉，所以需要通過研究范圍的邊界線把非封閉的等水位線封閉起來，以獲得封閉的等水位線。

②封閉的等水位線：在對沒有封閉的等水位面進行封閉后，等值面是通過該線與線之間的拓撲關系構建的。再對構建的等值面實施包含分析，如果面C與其他面A、面B存在包含關系，則將面A與面B進行疊加，再去刪減C；如果面C和其他面不存在包含關系，則無需刪減。如圖6和圖7所示。

Figure 6 Having an inclusion relationship圖6 存在包含關系

Figure 7 Not having an inclusion relationship圖7 不存在包含關系

以此類推，遍歷研究范圍內所有等值面，最后統計匯總經過處理后的面的大小，將面積不為零的區域用顏色渲染。顏色序列由淺到深表示水位埋深越來越大。

本文選取江蘇省鹽城市2013年7月的第Ⅲ承壓含水層水位監測數據，并繪制研究范圍等水位面與流線，結果如圖8和圖9所示。可以通過等水位面圖得到：在全局上，該市東部地區的地下水位明顯要比該市西部地區偏高，此外，該市市區的地下水水位普遍較低。圖10將向量線式流場與水位等值面圖進行疊加，可以通過流線的箭頭方向，更直觀地看出該市地下水的流向及其地下水供給關系。對水位紅線分析，可以清楚地看出哪一塊是超限采區，哪一塊是超禁采區，灰色區域為水位低于禁采水位紅線的范圍，較深灰色區域即為水位低于限采水位紅線的范圍。進一步分析水位紅線分析的結果，統計低于水位紅線的區域，可明確得知哪些地區的水位低于紅線，需要控制開采量。

Figure 8 Water level contour map圖8 水位等值面圖

Figure 9 Arrowhead representation of flow field圖9 點箭頭式表示流場

Figure 10 Superposition of vector line flow field and water level isogram圖10 向量線式流場與水位等值圖疊加

5 結束語

本文采用GIS技術設計的基于二維流場的地下水位紅線分析算法，分別從4個方面實現地下水流場的可視化，設計了一種新的種子點布局方法來表征地下水流場，采用改進的EULER方法進行數值積分來表達流線，將地下水流場和水位紅線管理進行綜合。討論了水位和紅線的幾何關系，總結了基于等水位面的生成方法，對研究范圍內地下水位空間分布特征實施科學可視化管理，利用線和面的空間拓撲聯系，將結論定量化分析。目前，等水位線和等水位面生成方法已有很多的科研成果，有效保證了水位紅線分析結果的精度。利用GIS符號化功能可使紅線分析結果更加直觀明確，將紅線分析結果與其他數據進行疊加，以得出更多有效結論，比如將紅線分析結果和研究范圍的行政區邊界、具體開采井點等數據進行疊加分析，能夠準確地觀察各地區的地下水開采狀況，并為地下水開采方案和地下水恢復管理對策的制定提供科學有效的數據支撐。