陳曉樂，崔益安，謝 靜，陽 兵，柳建新

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長沙 410083；2.中南大學(xué) 有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長沙 410083；3.中南大學(xué) 教育部有色金屬成礦預(yù)測重點(diǎn)實(shí)驗(yàn)室，長沙 410083)

0 引言

自然電位法是一種傳統(tǒng)的地球物理勘探方法，它利用地表或井中自然電位異常識別地質(zhì)與地球物理目標(biāo)體，在礦產(chǎn)、油氣勘探及垃圾填埋檢測等領(lǐng)域中都有成功應(yīng)用。自然電位法作為一種被動(dòng)源法，其野外觀測快捷，成本低廉。同時(shí)該方法對多種與污染擴(kuò)散伴生的滲流運(yùn)動(dòng)、氧化還原反應(yīng)以及微生物活動(dòng)信號較為敏感[1-3]，適用于工程與環(huán)境物探領(lǐng)域，尤其是地下有機(jī)污染物檢測及監(jiān)測。Naudet[1-2]和Revil[3]通過自然電位法圈定污染位點(diǎn)的氧化還原前沿(redox fronts);Arora[4]和Gallas[5]使用該方法成功監(jiān)測了垃圾填埋場的泄漏;Mao[6]通過沙箱和現(xiàn)場試驗(yàn)研究了自然電位層析成像，并將其應(yīng)用于地下水修復(fù)和污染羽流圈定;Abbas[7]分析了通過自然電位數(shù)據(jù)反演得到的一個(gè)富有機(jī)物污染場地的氧化還原場分布;Cui[8-9]進(jìn)行了沙箱實(shí)驗(yàn)，嘗試使用自然電位法來監(jiān)測地下污染的擴(kuò)散和遷移。

有機(jī)污染區(qū)域觀測到的地表自然電位異常跟地電池(Geobattery)模型的自然電位很相似[1-2，10]，但直接用地電池模型難以解釋這種自然電位異常，原因在于像垃圾填埋場這類有機(jī)污染區(qū)域不存在地電池模型中的電子導(dǎo)體礦物來傳遞電子。因此有研究人員提出可能存在微生物提供電子傳導(dǎo)通道并參與形成地電池模型的情況[1-4，10]。在此猜想下，Revil[3]基于地電池模型提出了生物地電池(Biogeobattery)模型來解釋有機(jī)污染區(qū)域的自然電位異常,此后生物地電池模型被廣泛應(yīng)用;Fachin[11]開展了物理模擬實(shí)驗(yàn)，在實(shí)驗(yàn)沙箱中埋設(shè)有機(jī)污染物成功地形成了生物地電池;Doherty[12]采用生物地電池模型來解釋在一個(gè)廢棄煤氣廠的煤焦油有機(jī)污染物引起的自然電位異常;Risgaard-Petersen[13]還將生物地電池模型用于海洋淤泥和海底沉積物中的生物電化學(xué)過程研究。以上研究表明，地電池多源模型能有效表征地下復(fù)雜氧化還原反應(yīng)的電位異常特性。基于地電池多源模型，開展地下污染物氧化還原場數(shù)值模擬研究，有助于自然電位觀測數(shù)據(jù)的反演解釋，更好地發(fā)揮自然電位法在環(huán)境領(lǐng)域的應(yīng)用效果。

無限單元用于解決截?cái)噙吔鐔栴}，于上世紀(jì)70年代提出[14]。無限單元可分為三類：①Astley映射無限單元[15]；②Bettess映射無限單元[16]；③Burnett多極開展無限單元[17]。與有限單元相比，無限單元更合適求解無限域或半無限域問題。肖曉[18]、原源[19]實(shí)現(xiàn)了二維有限元-無限元耦合法，并將其成功應(yīng)用于2.5維直流電阻率正反演研究中;湯井田[20]、張林成[21]、Xie[22-23]先后開展了三維有限元-無限元耦合法研究，并將其應(yīng)用于地球物理數(shù)值模擬中。常規(guī)有限元施加截?cái)噙吔鐥l件比較麻煩，而且不利于求解場源靠近邊界的情況。早期二維有限元-無限元耦合法的實(shí)現(xiàn)均基于矩形格[18-19]，對于復(fù)雜模型的剖分效果不佳，適應(yīng)性不強(qiáng)。這里基于三角形有限元及應(yīng)用最為廣泛的Astley型單向映射無限元，改進(jìn)了二維有限元-無限元耦合算法，通過均勻半空間點(diǎn)電源模型驗(yàn)證該算法的正確性及有效性，并將其應(yīng)用于地下有機(jī)污染物降解過程所產(chǎn)生的氧化還原場的數(shù)值模擬研究中，探討自然電位法在環(huán)境污染探查、監(jiān)測領(lǐng)域的可行性。數(shù)值結(jié)果表明，筆者所改進(jìn)的算法行之有效，而自然電位法在環(huán)境污染領(lǐng)域也有著其獨(dú)特的勘探優(yōu)勢。

1 2.5維自然電位邊值問題

首先考慮二維模型的邊值問題，在二維地電模型中，電流源仍為三維源，故需利用空間波數(shù)域的狄拉克函數(shù)特性求解2.5維自然電位問題。2.5維自然電位地電模型滿足的邊值問題為[18-19]：

▽·(σ▽U)-k2σU=-I0δ(A) ∈Ω

(1)

(2)

(3)

其中：σ為介質(zhì)電導(dǎo)率；I0為點(diǎn)源電流強(qiáng)度；δ(A)為與點(diǎn)源位置A有關(guān)的δ函數(shù)；n為邊界外法向單位矢量；k為離散波數(shù)；U為波數(shù)域電位；K0(kr)、K1(kr)分別為零階、一階第二類修正貝塞爾函數(shù)；Ω、Γs、?！薹謩e為經(jīng)傅里葉變換后的研究區(qū)域、地表邊界、無窮遠(yuǎn)邊界；r為點(diǎn)源點(diǎn)到外邊界任意點(diǎn)的距離(圖1)。

圖1 參數(shù)示意圖Fig.1 The sketch of parameters

在有限元-無限元耦合法中，無需考慮式(3)，而由無限單元充當(dāng)邊界條件。故有限元-無限元耦合法的邊值問題所對應(yīng)的變分問題為：

(4)

δF(U)=0

(5)

式(4)右端的變分為：

(6)

2 有限元-無限元耦合法

有限元-無限元的空間耦合如圖2所示。在有限元區(qū)域，利用三角形有限元的3個(gè)節(jié)點(diǎn)構(gòu)造有限單元形函數(shù)，采用該組形函數(shù)φi近似計(jì)算點(diǎn)x=(x,z)的波數(shù)域電位為式(7)。

圖2 有限元-無限元的空間耦合Fig.2 The spacial coupling of finite-infinite elements

(7)

將式(7)代入式(6)有：

(8)

對式(8)中的面積分項(xiàng)做變換有：

(9)

對于無限元區(qū)域，采用的是Astley型單向映射無限單元對其進(jìn)行剖分[15，18-19]。該無限單元的具體映射過程如圖3所示。

圖3 二維四節(jié)點(diǎn)無限元映射過程Fig.3 The mapping process of 2D infinite elements(a)母單元;(b)子單元

在二維四節(jié)點(diǎn)無限單元中，任意位置的空間坐標(biāo)可表示為：

(10)

(11)

其中:Mi為坐標(biāo)映射函數(shù)，其表達(dá)式為：

(12)

(13)

(14)

(15)

式中：ξi，ηi為局部坐標(biāo)。經(jīng)過上述坐標(biāo)映射之后，電位借助無限單元延伸至無限遠(yuǎn)處，并且在無限遠(yuǎn)處衰減為零。

無限單元中任意位置的波數(shù)域電位可表示為：

(16)

其中：Ui為各節(jié)點(diǎn)電位；Ni為無限元插值形函數(shù)，其表達(dá)式為：

(17)

(18)

(19)

(20)

無限元區(qū)域所滿足的偏微分方程及其變分形式，同樣可分別由式(1)、式(2)、式(4)、式(5)和式(8)表示。故無限元區(qū)域的面積分可表示為式(21)。

(21)

式(9)、式(21)中的Φ分別由有限單元及無限元求取。合并式(9)、式(21)，有：

δUT[(K1+K2)U-F]=0

(22)

其中

(23)

(24)

(25)

式(23)與式(24)的耦合過程可表示為[22-23]：

(26)

3 傅里葉正反變換

這里所采用的波數(shù)及其權(quán)值如表1所示。

由式(22)、式(26)有：

KU=F

(27)

其中:K=K1+K2;U為波數(shù)域電位向量;F為源向量。

在不同波速條件下(表1)，求解方程(27)得到波數(shù)域電位。然后由傅里葉反變換求解得到空間域電位：

(28)

其中,u為空間域電位。經(jīng)以上計(jì)算，從而得到二維模型中所有節(jié)點(diǎn)的電位分布。

4 數(shù)值算例

4.1 算法可靠性驗(yàn)證

采用均勻半空間多源模型進(jìn)行算法可靠性驗(yàn)證，模型如圖4所示。模型尺度為40 m×40 m，剖分為40×40個(gè)四邊形網(wǎng)格。為具有更好的對比性，在有限元-無限元耦合法中，通過將四邊形網(wǎng)格對半剖分，實(shí)現(xiàn)簡單的三角形剖分。點(diǎn)電流源幅值均為1 A，波數(shù)及其權(quán)值選擇如表1所示。數(shù)值模擬結(jié)果(圖5)表明，在相同的網(wǎng)格剖分(節(jié)點(diǎn)分布)條件下，有限元-無限元耦合法的計(jì)算精度與具有混合邊界條件的常規(guī)有限單元法相當(dāng)，相對于解析解的均方誤差分別為0.022 7、0.019 0。而實(shí)施過程中無需計(jì)算邊界積分，計(jì)算效率更高，尤其在多源動(dòng)態(tài)模型中更能進(jìn)一步體現(xiàn)其優(yōu)勢，驗(yàn)證了算法的正確性及有效性。

表1 模型計(jì)算所采用的波數(shù)k及其權(quán)值g[24]Tab.1 Wave number k and its weight g used in the model calculations

圖4 均勻半空間多源模型Fig.4 Homogeneous half-space multi-source model

圖5 地表電位曲線Fig.5 Surface electric potential curves

4.2 氧化還原地電模型

地下有機(jī)污染物的降解可視為動(dòng)態(tài)過程[25-26]。相關(guān)生物電化學(xué)反應(yīng)所構(gòu)建的地電池模型是隨著時(shí)間動(dòng)態(tài)變化的。對于每個(gè)靜態(tài)時(shí)刻，其地電模型可由圖7近似表征。隨著氧化還原強(qiáng)度的變化，圖7模型中的電流強(qiáng)度發(fā)生動(dòng)態(tài)變化。在實(shí)際問題中，地下有機(jī)污染物的降解進(jìn)程受到許多因素的影響，如土壤中自由空氣的濃度、有機(jī)污染物的種類及其濃度、土壤微生物的類型及其數(shù)量分布等。但結(jié)合微生物的周期性代謝特征，整個(gè)降解系統(tǒng)的氧化還原反應(yīng)強(qiáng)度的總體趨勢可簡化為：先緩后急，再緩至峰值，最后再逐漸衰減。根據(jù)該規(guī)律，筆者由圖6所示的假設(shè)曲線近似表征整個(gè)降解過程的電流源幅值變化，以定性探究該過程的自然電位異常特征。

地電模型(圖7)可近似模擬微生物降解地下有機(jī)污染物時(shí)在地下水位附近發(fā)生的氧化還原反應(yīng)所引起的自然電位異常。模型中有限單元區(qū)域尺度為20 m×30 m，采用三角形剖分，除地表外的三個(gè)邊界由無限單元作為邊界單元充當(dāng)邊界條件，為四邊形剖分。設(shè)定地下水位位于2 m深度處，上層電阻率為10 Ω·m，下層電阻率為100 Ω·m。假定地下水位附近有兩處位置存在微生物降解有機(jī)污染物的現(xiàn)象，其中地下水位上部發(fā)生氧化反應(yīng)，積累負(fù)電荷，地下水位下部發(fā)生還原反應(yīng)，積累正電荷，從而形成生物地電池模型。通過離散正負(fù)點(diǎn)電流源近似正負(fù)電荷的區(qū)域性分布，各點(diǎn)電流源幅值的動(dòng)態(tài)變化如圖6所示。

圖6 電流強(qiáng)度幅值的動(dòng)態(tài)變化Fig.6 The dynamic variation of current intensity amplitudes

圖7 生物地電池簡化模型Fig.7 The simplified model of biogeobattery

數(shù)值模擬結(jié)果(圖8、圖9)表明，地下有機(jī)污染物降解過程能在地面觀測到自然電位負(fù)異常，且通過局部負(fù)異常極值點(diǎn)能判斷氧化還原反應(yīng)的分布范圍及場點(diǎn)數(shù)量。地表自然電位曲線也能進(jìn)一步反映地下有機(jī)污染物降解過程的變化情況，與圖6的變化趨勢呈正相關(guān)。

圖8 異常及非異常點(diǎn)自然電位隨電流源的變化Fig.8 The variation of self-potential of abnormal and non-abnormal points with current sources

圖9 地表自然電位異常曲線Fig.9 Surface self-potential curves

5 結(jié)論

改進(jìn)了原有基于矩形網(wǎng)格的二維有限元-無限元耦合算法，將有限元區(qū)域三角化，更適合模擬復(fù)雜地電模型。首先通過均勻半空間模型驗(yàn)證了該算法的可靠性，然后將其應(yīng)用于微生物降解地下有機(jī)污染物過程所產(chǎn)生的氧化還原場的數(shù)值模擬中，研究地表自然電位對該過程的異常響應(yīng)。得到以下結(jié)論：

1) 基于有限元區(qū)域三角剖分的改進(jìn)算法能更有效地剖分復(fù)雜地質(zhì)模型，在相同的網(wǎng)格剖分條件下，單次計(jì)算速度及計(jì)算精度與常規(guī)有限單元法相當(dāng)，且方便處理場源信息及邊界條件，更適用于自然電位等多源動(dòng)態(tài)模型的數(shù)值模擬。

2) 地質(zhì)微生物參與的地下有機(jī)污染物降解過程會在地表呈現(xiàn)自然電位負(fù)異常，且可通過地表電位的變化情況評估地下有機(jī)污染物的降解進(jìn)度。

3) 動(dòng)態(tài)模型的模擬結(jié)果表明地表自然電位異常與地下電流強(qiáng)度呈正相關(guān)。