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

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

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

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

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

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

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

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

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

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

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

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

圖2 六面體單元剖分為5個(gè)四面體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 三維點(diǎn)源場邊值問題

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

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

(1)

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

在地表Г上滿足

?u/?n=0

(2)

在近似無窮遠(yuǎn)邊界Г∞上，取混合邊界條件

(3)

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

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

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

2.2 四面體單元分析

與三維點(diǎn)源模型邊值問題等價(jià)的變分問題為

(4)

在四面體單元中(圖7)，空間坐標(biāo)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具體表達(dá)式見文獻(xiàn)[14]，k=i,j,l,m；而Ve為四面體體積。

對(duì)全域進(jìn)行剖分，將變分問題的積分分解為各單元e上的積分。對(duì)體積分項(xiàng)有

(6)

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

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

(7)

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

對(duì)全體節(jié)點(diǎn)組成的矩陣有

(8)

3 正演模擬

3.1 算法驗(yàn)證

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

3.2 誤差分析

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

3.3 滲流自然電場地電模型

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

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

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

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

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

傾斜滲流模型如圖13所示。模型大小為100 m×100 m ×100 m，滲流口位于模型中部，以3個(gè)不同傾角(0°、30°、45°)向地下滲流。以正負(fù)離散點(diǎn)源近似模擬滲流場的電荷分布，模擬結(jié)果如圖14、圖15。傾斜滲流走向會(huì)引起地表電位的不對(duì)稱分布：在滲流方向上能觀測到自然電位正異常，且隨著傾角的增大，即滲流通道越趨于水平，所觀測到的正異常越大。據(jù)此可有效判斷地下滲流走向，可應(yīng)用于地下水監(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)缺點(diǎn)，并提出一種新的三維結(jié)構(gòu)化網(wǎng)格剖分技術(shù)；然后結(jié)合三維點(diǎn)源場邊值問題，推導(dǎo)了四面體有限單元基本方程，并對(duì)比分析了幾種剖分技術(shù)的模擬精度；最后實(shí)現(xiàn)了滲流自然電場地電模型正演模擬，探討了滲流走向以及正負(fù)地形對(duì)地表自然電位分布的影響。

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

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