自然電場三維有限元正演模擬

周竹生,朱海倫,謝 靜,劉思琴,楊 陽

(1．中南大學 地球科學與信息物理學院,長沙 410083；2．山西省煤炭地質(zhì)物探測繪院,山西 晉中 030600)

起伏地形對地表電位分布的影響較大，研究起伏地形引起的電位假異常是電法勘探理論發(fā)展的需要[1-4]。目前直流電法二維正反演技術(shù)已相當成熟，也開展了大量的三維正反演研究工作并取得了顯著成果；但復雜地形條件下的三維正反演研究仍有探索空間。正演是反演的基礎，開展快速、高精度的正演模擬工作，有助于反演技術(shù)的發(fā)展。常規(guī)的正演模擬方法，如有限單元法、有限體積法、有限差分法等，通過簡化實際地電模型實現(xiàn)正演模擬過程，以更好地認識地下異常體的地表異常響應。其中有限單元方法模擬精度高、求解簡易且規(guī)范、自動滿足內(nèi)部邊界條件，適于各種復雜的物性分布情況，是一種實用高效并應用廣泛的正演模擬方法[5-9]。

單元網(wǎng)格剖分是正演模擬的核心工作，研究網(wǎng)格剖分方法對建立有限元模型至關(guān)重要。目前常規(guī)的結(jié)構(gòu)化網(wǎng)格剖分技術(shù)仍以六面體為基礎單元。規(guī)則六面體剖分只適用于水平地形條件下的三維正演模擬，一定程度上限制了三維有限元的發(fā)展和應用。熊彬等首先提出先將研究區(qū)域進行一級六面體剖分，再將六面體二級剖分為6個四面體，對四面體添加地形數(shù)據(jù)，以完成三維復雜地形條件下的有限元正演模擬[10]；在此基礎上，呂玉增等提出一種四面體網(wǎng)格交叉剖分技術(shù)，并將六面體的二級剖分減少為5個四面體[11]。但以上兩種剖分技術(shù)必須滿足四面體單元交叉分布，才能保證模擬誤差呈對稱性分布。強建科提出任意三棱柱單元剖分方式，該三棱柱包含地形特征，能有效貼合起伏地形，成功實現(xiàn)復雜地形條件下的三維有限元正演模擬[12]；但該剖分技術(shù)必須滿足模型的頂?shù)捉缑嫫叫?，不利于模型的建立，且對模擬精度有一定影響。

本文基于前人研究，提出一種新的三維有限元結(jié)構(gòu)化網(wǎng)格剖分技術(shù)，該技術(shù)先將研究區(qū)域一級剖分為規(guī)則六面體，再將六面體二級剖分為4個三棱柱，最后將每個三棱柱三級剖分為3個四面體。該剖分技術(shù)綜合了三棱柱剖分和四面體剖分的優(yōu)點，能更好地貼合復雜地形，且有效解決了剖分方式引起的模型底界面畸變以及誤差不對稱分布等問題；同時在研究區(qū)域一定、六面體網(wǎng)格尺度一定時，模擬得到的節(jié)點信息更為豐富，模擬精度更高。

本文首先通過均勻半空間三維點源模型進行算法驗證，同時建立滲流自然電場地電模型，探討分析滲流走向以及起伏地形對地表自然電位分布的影響。

1 網(wǎng)格剖分技術(shù)

1.1 常規(guī)結(jié)構(gòu)化網(wǎng)格剖分技術(shù)

常規(guī)六面體網(wǎng)格剖分技術(shù)是最簡單易行的三維剖分方式，但其最大缺陷是不能實現(xiàn)復雜地形條件下的三維模擬。熊彬等[10]將研究區(qū)域先剖分為規(guī)則六面體，然后將每個六面體單元再二次剖分為6個四面體單元以添加高程數(shù)據(jù)(圖1)；呂玉增等[11]提出的四面體網(wǎng)格交叉剖分技術(shù)是將六面體網(wǎng)格中的四面體網(wǎng)格減少到5個(圖2)；強建科[12]提出的任意三棱柱網(wǎng)格剖分技術(shù)能較好地貼合起伏地形，但每個三棱柱單元都必須滿足上、下底面平行相等(圖3)，故而該剖分技術(shù)會致使整個研究區(qū)域的頂?shù)走吔缤耆叫邢嗟?，從而使得整個研究區(qū)域的空間結(jié)構(gòu)變得畸形，將對模擬精度產(chǎn)生一定影響。以上剖分技術(shù)有一個共同的缺陷，即每個六面體單元的上頂面劃分為2個三角形時(圖4)，會導致計算誤差分布不具對稱性。雖然呂玉增等[11]提出的交叉剖分技術(shù)能解決對稱精度問題，但該過程不易于程序?qū)崿F(xiàn)。

1.2 新結(jié)構(gòu)化網(wǎng)格剖分技術(shù)

綜合考慮以上幾種常規(guī)的三維有限元結(jié)構(gòu)化網(wǎng)格剖分技術(shù)的優(yōu)缺點，本文提出一種新的剖分技術(shù)，集以上所有剖分方式的優(yōu)點于一身，同時摒棄了上述剖分方法的所有缺陷。該剖分方式是將1個六面體先剖分為4個三棱柱，然后再將各三棱柱二次剖分為3個四面體，即將1個六面體剖分為12個四面體，相應單元節(jié)點數(shù)由8增加到10(圖5)。相比傳統(tǒng)的結(jié)構(gòu)化網(wǎng)格剖分技術(shù)，新剖分技術(shù)在六面體頂?shù)酌娓髟黾恿艘粋€中心節(jié)點，在研究區(qū)域一定、六面體剖分網(wǎng)格一定時，能有效提高正演模擬精度，同時保證了計算誤差呈對稱性分布，且不用考慮復雜的交叉剖分問題，便于編程實現(xiàn)。

圖1 六面體單元剖分為6個四面體Fig.1 A hexahedral element is subdivided into six tetrahedrons

圖2 六面體單元剖分為5個四面體Fig.2 A hexahedral element is subdivided into five tetrahedrons

圖3 六面體單元剖分為三棱柱網(wǎng)格單元Fig.3 Hexahedral elements are subdivided into triangular prism elements

圖4 常規(guī)六面體上頂面剖分方式Fig.4 Top surface subdivision of hexahedron elements

2 有限單元法

2.1 三維點源場邊值問題

均勻半空間三維點源模型(圖6)，電位滿足微分方程[13]

▽·(σ▽u)=-2Iδ(A)

(1)

其中：σ為電導率；u為電位；I為點源電流；δ(A)為關(guān)于場源所在位置A點的δ函數(shù)。

在地表Г上滿足

?u/?n=0

(2)

在近似無窮遠邊界Г∞上，取混合邊界條件

(3)

其中：n為邊界外法向單位向量；r為點源指向邊界Г∞上任一點的向量，而區(qū)域內(nèi)部電導率邊界為自然邊界條件，不予考慮。

圖5 新剖分技術(shù)的剖分過程Fig.5 The subdivision process of the new subdivision technology

圖6 均勻半空間點源模型Fig.6 3D point source model in homogeneous half-space

2.2 四面體單元分析

與三維點源模型邊值問題等價的變分問題為

(4)

在四面體單元中(圖7)，空間坐標x、y、z及電位u的插值函數(shù)為

圖7 四面體單元Fig.7 Tetrahedral element

(5)

式中：Ψ分別為x、y、z；Nk為形函數(shù)，且Nk=(ak+bkx+cky+dkz)/(6Ve)，其中ak、bk、ck、dk具體表達式見文獻[14]，k=i,j,l,m；而Ve為四面體體積。

對全域進行剖分，將變分問題的積分分解為各單元e上的積分。對體積分項有

(6)

整理得到四面體單元體積分矩陣為：K1e=(kij)，其中kij=σ(bibj+cicj+didj)/(36Ve)。

對于邊界積分，近似cos(r,n)及|r|為常量，將其提到積分號外，整理得到邊界積分為

(7)

其中：K2e為邊界積分系數(shù)矩陣[14]；ue為單元節(jié)點電位向量。

對全體節(jié)點組成的矩陣有

(8)

3 正演模擬

3.1 算法驗證

本節(jié)通過建立均勻半空間點源模型進行算法驗證。模型空間剖分為30×30×30個六面體網(wǎng)格單元，區(qū)域大小取60 m×60 m×30 m，點源置于地表中心。將各個六面體單元分別按圖4的2種剖分方式以及本文提出的新剖分技術(shù)進行剖分，實現(xiàn)不同剖分方式的正演模擬?？紤]點源奇異性，圖8所示為地下2 m深處的電位等值線平面，正演結(jié)果表明剖分方式對計算精度有一定影響，單一朝向的剖分技術(shù)會使得本應呈對稱性分布的地表電位失去對稱性。

3.2 誤差分析

為進一步說明新剖分技術(shù)對模擬精度的影響，將圖8中的3個電位平面分別與解析解作誤差分析，剖分方式的不同導致誤差分布存在較大差異(圖9)。單一朝向的剖分方式使得誤差不再關(guān)于點源呈對稱分布，而是具有方向性(呈軸對稱分布)。本文提出的剖分方式，能保證誤差對于點源呈對稱分布，且模擬精度更高。

3.3 滲流自然電場地電模型

自然電場法無需供電，儀器設備輕便，生產(chǎn)效率較高，在傳統(tǒng)的礦產(chǎn)與油氣勘探[15]、地熱火山調(diào)查[16-17]、地震觀測[18-19]等領域得到廣泛應用，近年來還在工程與環(huán)境[20-24]等新興領域得到廣泛應用。流動電位是基于電動機理產(chǎn)生的自然電位，其在自然條件下普遍存在，在工程和地下水地球物理學中受到廣泛關(guān)注。它為檢測水壩、水庫地面、排水結(jié)構(gòu)滲漏等問題提供了一種經(jīng)濟有效的手段，還可以用于監(jiān)測含水層中的水運動、山體滑坡、垃圾填埋場污染泄漏等問題。有針對性地開展流動電位正演模擬研究將有助于反演解釋，更好地提高該方法的應用效果。利用本文提出的結(jié)構(gòu)化網(wǎng)格剖分技術(shù)，實現(xiàn)流動電位三維正演模擬，探討分析起伏地形以及滲流方向引起的地表自然電位異常特征。

圖8 深度為2 m處的電位等值線圖Fig.8 Electric potential contour maps at the depth of 2 m(A)六面體頂面向左傾斜;(B)六面體頂面向右傾斜;(C)新剖分技術(shù)

圖9 不同剖分方式的相對誤差Fig.9 Relative error of different subdivision techniques(A)六面體頂面向左傾斜;(B)六面體頂面向右傾斜;(C)新剖分技術(shù)

直立滲流模型如圖10所示。2個模型大小均為100 m×100 m×100 m，滲流通道位于中部，在模型(B)中添加正負地形。通過正負離散點源近似模擬滲流場的電荷分布，數(shù)值模擬結(jié)果如圖11、圖12。正負地形對地表電位分布均有一定影響：因正地形使得電位衰減距離變長，故異常幅度減?。欢摰匦蔚刃в陔娢凰p距離變短，故異常振幅增大。

圖10 三維模型的二維剖面Fig.10 2D subsurface structures of 3D models(A)直立滲流無地形；(B)直立滲流帶正負地形

傾斜滲流模型如圖13所示。模型大小為100 m×100 m ×100 m，滲流口位于模型中部，以3個不同傾角(0°、30°、45°)向地下滲流。以正負離散點源近似模擬滲流場的電荷分布，模擬結(jié)果如圖14、圖15。傾斜滲流走向會引起地表電位的不對稱分布：在滲流方向上能觀測到自然電位正異常，且隨著傾角的增大，即滲流通道越趨于水平，所觀測到的正異常越大。據(jù)此可有效判斷地下滲流走向，可應用于地下水監(jiān)測、追蹤地下污染滲流等。

圖11 地表自然電位等值線圖Fig.11 Contour maps of ground self-potential

圖12 主剖面地表自然電位曲線Fig.12 Ground self-potential curves of main profile

圖13 傾斜滲流模型的二維結(jié)構(gòu)Fig.13 2D subsurface structure of inclined seepage model

圖14 地表自然電位等值線圖Fig.14 Contour maps of ground self-potential (A)傾角0°；(B)傾角30°；(C)傾角45°

圖15 主剖面地表電位曲線Fig.15 Ground self-potential curves of main profiles

4 結(jié) 論

本文從網(wǎng)格剖分技術(shù)出發(fā)，分析總結(jié)了幾種常用的三維結(jié)構(gòu)化網(wǎng)格剖分技術(shù)的優(yōu)缺點，并提出一種新的三維結(jié)構(gòu)化網(wǎng)格剖分技術(shù)；然后結(jié)合三維點源場邊值問題，推導了四面體有限單元基本方程，并對比分析了幾種剖分技術(shù)的模擬精度；最后實現(xiàn)了滲流自然電場地電模型正演模擬，探討了滲流走向以及正負地形對地表自然電位分布的影響。

在算法方面，該剖分技術(shù)具有網(wǎng)格單元呈對稱性分布的特點，同時易編程實現(xiàn)，不用考慮網(wǎng)格交叉技術(shù)的繁瑣過程；在正演模擬結(jié)果方面，該剖分技術(shù)能完美解決常規(guī)結(jié)構(gòu)化網(wǎng)格剖分技術(shù)中的誤差非對稱性問題，同時在一定程度上提高了模擬精度。由滲流自然電場正演模擬結(jié)果可知，正地形使電位異常幅度減小，而負地形使電位異常振幅增大；在滲流方向上能觀測到較大的正異常，有助于確定滲流走向。

需要說明的是，本文提出的剖分技術(shù)在一定程度上增加了單元節(jié)點數(shù)，增加了計算量，但計算精度也更高。當今計算機性能高速發(fā)展，以較少的計算增加量換取更高的計算精度是有一定實際意義的，在選擇剖分技術(shù)時，綜合考量側(cè)重點即可。