基于Debye模型的時(shí)間域航空電磁激電效應(yīng)正演模擬

張?chǎng)纬?殷長(zhǎng)春

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130026

0 引言

時(shí)間域航空電磁法采用機(jī)載平臺(tái)搭載電磁儀器,可以在沙漠、高山、沼澤、湖泊等地形條件復(fù)雜的地區(qū)進(jìn)行勘探,是一種高效的地球物理手段,現(xiàn)已被廣泛應(yīng)用于礦產(chǎn)資源、地下水及環(huán)境工程等眾多勘探領(lǐng)域[1]。與頻率域電磁法相比,時(shí)間域航空電磁法對(duì)深部目標(biāo)具有更好的分辨率和探測(cè)靈敏度[2]。

隨著時(shí)間域航空電磁系統(tǒng)的進(jìn)步和數(shù)據(jù)采集精度的不斷提高,地球物理工作者們經(jīng)常會(huì)遇見晚期道時(shí)間域航空電磁響應(yīng)出現(xiàn)符號(hào)反轉(zhuǎn)的現(xiàn)象,這種現(xiàn)象多發(fā)生在地下含硫化物以及一些含水地區(qū)。學(xué)者們認(rèn)為該現(xiàn)象是由于激發(fā)極化效應(yīng)引起的[3-5]。Weidelt[6]也進(jìn)一步通過理論證明,對(duì)于共中心回線裝置,時(shí)間域電磁響應(yīng)數(shù)據(jù)晚期時(shí)間道中的符號(hào)反轉(zhuǎn)現(xiàn)象出現(xiàn)的原因往往是激發(fā)極化效應(yīng)。然而,傳統(tǒng)的僅依賴電阻率的解釋方法無法對(duì)該現(xiàn)象進(jìn)行有效反演,給時(shí)間域航空電磁數(shù)據(jù)解釋帶來了一定的困擾。因此,時(shí)間域航空電磁法激電效應(yīng)研究具有重要的意義。

近年,時(shí)間域航空電磁數(shù)據(jù)成像和一維反演解釋技術(shù)已經(jīng)趨于成熟。然而,實(shí)際情況下大地多為復(fù)雜的三維結(jié)構(gòu),一維成像和反演往往難以得到地下真實(shí)的電性分布結(jié)構(gòu)[7-8]。因此,發(fā)展三維反演成為必然趨勢(shì)。另一方面,正演是反演的基礎(chǔ),為對(duì)航空電磁數(shù)據(jù)進(jìn)行有效解釋,需要精確的正演算法。因此,本文將重點(diǎn)研究時(shí)間域航空電磁三維激電效應(yīng)正演模擬。

針對(duì)激電效應(yīng)的正演算法目前主要分為兩大類。一種是基于頻時(shí)轉(zhuǎn)換方法[9-11],另一種是對(duì)麥克斯韋方程在空間和時(shí)間域直接進(jìn)行離散[2, 12-13]。然而,傅里葉變換的精度受采樣和變換方法的影響很大[14],因此直接求解時(shí)間域麥克斯韋方程受到了更多的關(guān)注。

在激發(fā)極化介質(zhì)中,電導(dǎo)率與頻率有關(guān)。Cole兄弟[15]在電化學(xué)研究中引入了復(fù)介電常數(shù)的表達(dá)式來描述液體的頻散情況;隨后Pelton等[16]在前人研究基礎(chǔ)上,通過大量巖礦石測(cè)量結(jié)果,提出了國(guó)內(nèi)外學(xué)者廣泛采用的Cole-Cole模型,通過數(shù)學(xué)模型來描述均勻巖礦石由激電效應(yīng)引起的復(fù)電阻率隨頻率的變化。目前,除了Cole-Cole模型[17-19],還有Debye模型[20-22]和Cole-Davidson模型[23]。其中,Debye模型是Cole-Cole模型的特例(頻率相關(guān)系數(shù)c=1)[24],數(shù)學(xué)形式簡(jiǎn)單、最具代表性。本文選用Debye模型來描述巖礦石激電效應(yīng)引起的復(fù)電阻率隨頻率的變化,進(jìn)而對(duì)時(shí)間域航空電磁激電效應(yīng)進(jìn)行三維正演模擬。

在空間離散方面,本文采用基于非結(jié)構(gòu)四面體網(wǎng)格的矢量有限元方法。該方法可以有效克服復(fù)雜幾何形體和復(fù)雜起伏地形因建模造成的誤差[25-26]。在時(shí)間離散方面,本文選用無條件穩(wěn)定的隱式二階后推歐拉離散格式。首先從麥克斯韋方程出發(fā),采用二階后推歐拉離散格式對(duì)傳導(dǎo)電流和源電流密度時(shí)間導(dǎo)數(shù)進(jìn)行時(shí)間離散;然后將Debye模型引入到歐姆定律中,推導(dǎo)出時(shí)間域歐姆定律表達(dá)式,進(jìn)而利用二階后推歐拉離散格式對(duì)歐姆定律中的時(shí)間導(dǎo)數(shù)進(jìn)行離散,再利用非結(jié)構(gòu)有限元方法對(duì)空間進(jìn)行離散,建立方程進(jìn)行求解;最后通過典型模型驗(yàn)證方法的準(zhǔn)確性,并對(duì)異常響應(yīng)特征進(jìn)行分析。

1 方法原理

1.1 控制方程

時(shí)間域麥克斯韋方程可以寫成:

(1)

(2)

式中:e(r,t)和h(r,t)分別為r處t時(shí)刻的電場(chǎng)強(qiáng)度和磁場(chǎng)強(qiáng)度;μ為自由空間磁導(dǎo)率,本文假設(shè)為自由空間值4π×10-7H/m;j(r,t)為傳導(dǎo)電流密度;js(r,t)為源電流密度;d(r,t)為電位移矢量。式(2)中的右端第二項(xiàng)為位移電流項(xiàng)。航空電磁中由于位移電流遠(yuǎn)遠(yuǎn)小于傳導(dǎo)電流,因此可以忽略不計(jì)[27]。將式(1)(2)聯(lián)立并消去h,可得到如下方程:

(3)

對(duì)于式(3)中的時(shí)間導(dǎo)數(shù),我們采用無條件穩(wěn)定的二階后推歐拉離散格式進(jìn)行近似,則傳導(dǎo)電流密度和源電流密度的時(shí)間導(dǎo)數(shù)可分別離散為:

(4)

(5)

式中,時(shí)間步長(zhǎng)Δtk-1=tk-tk-1。則式(3)可整理為

(6)

1.2 極化介質(zhì)歐姆定律

傳導(dǎo)電流密度可由歐姆定律描述,即

j(r,t)=σe(r,t)。

(7)

式中,σ為電導(dǎo)率,在非極化介質(zhì)中是不變的。然而,當(dāng)存在激發(fā)極化效應(yīng)時(shí),電導(dǎo)率是頻散的,可用Debye模型[24]來描述:

(8)

式中:σ0為零頻電導(dǎo)率;m為充電率;τ為時(shí)間常數(shù);i為虛數(shù)單位;ω為角頻率。此時(shí),電流密度和電場(chǎng)強(qiáng)度不再滿足式(7)簡(jiǎn)單的線性關(guān)系。下面推導(dǎo)極化介質(zhì)中歐姆定律的形式。

首先我們給出極化介質(zhì)中頻率域歐姆定律,即

j(ω)=σ(ω)e(ω)。

(9)

則將式(8)代入到式(9)并整理可得

σ0e(ω)+iωτσ0e(ω)=j(ω)+iωτ(1-m)j(ω)。

(10)

應(yīng)用傅里葉逆變換,式(10)可變換到時(shí)間域[2],即

(11)

式(11)即為極化介質(zhì)中的歐姆定律,其左右兩端均出現(xiàn)了時(shí)間導(dǎo)數(shù)。依然采用二階后推歐拉離散格式對(duì)其進(jìn)行離散,則式(11)可整理為

(12)

為方便后續(xù)推導(dǎo),將式(12)相關(guān)系數(shù)用一些字符替代并整理為

(13)

其中:

將式(13)代入到式(6)中,即可得到電場(chǎng)的擴(kuò)散方程為

(14)

1.3 有限元方法

本文使用矢量有限元法對(duì)式(14)進(jìn)行空間離散。對(duì)于非結(jié)構(gòu)四面體網(wǎng)格,第n個(gè)網(wǎng)格r處tk時(shí)刻的電場(chǎng)強(qiáng)度可表示為

(15)

(16)

式中,j=1, 2, ???, 6。令加權(quán)殘差Rn=0,并通過一系列代數(shù)運(yùn)算,可以建立第n個(gè)網(wǎng)格的如下有限元方程:

(17)

式中,相關(guān)的矩陣可寫為:

(18)

將全部網(wǎng)格的式(18)組裝起來,可以得到tk時(shí)刻電場(chǎng)強(qiáng)度的稀疏對(duì)稱方程組,即

[2Δtk-1S+3g1M]e(tk)=12g2Me(tk-1)-3g2Me(tk-2)-12g3Mj(tk-1)+3g3Mj(tk-2)+4Mj(tk-1)-Mj(tk-2)-2Δtk-1Js。

(19)

通過對(duì)式(19)施加Dirichlet邊界條件n×e(tk)|boundary=0,可以得到最終的正演模擬方程為

Ae(tk)=b(tk)。

(20)

式中:A為系數(shù)矩陣;e(tk)為模型剖分網(wǎng)格單元各棱邊上tk時(shí)刻的待求電場(chǎng)強(qiáng)度;b(tk)為tk時(shí)刻的右端源項(xiàng)。

我們使用直接求解器MUMPS求解得到各棱邊的電場(chǎng)[28],而電場(chǎng)和磁場(chǎng)對(duì)時(shí)間的導(dǎo)數(shù)可以分別用矢量基函數(shù)及其導(dǎo)數(shù)通過插值計(jì)算得到。

2 數(shù)值算例

2.1 精度驗(yàn)證

為了驗(yàn)證本文算法的準(zhǔn)確性,我們計(jì)算了均勻半空間條件下不同激發(fā)極化參數(shù)組合模型的電磁響應(yīng),并將結(jié)果與一維半解析解進(jìn)行對(duì)比。本文采用中心回線裝置(圖1),系統(tǒng)參數(shù)設(shè)定如下:發(fā)射線圈(T)采用正八邊形,半對(duì)角線長(zhǎng)為15 m,發(fā)射線圈和接收線圈(R)高度均為30 m,發(fā)射電流為單位下階躍波。

圖1 航空電磁均勻半空間模型

空氣電導(dǎo)率設(shè)為10-8S/m,分別計(jì)算不同均勻半空間電導(dǎo)率、充電率、時(shí)間常數(shù)的電磁響應(yīng),并與一維半解析解進(jìn)行對(duì)比(圖2—圖4)。

a. 電磁響應(yīng);b. 與一維半解析解的相對(duì)誤差。σ0=0.010 S/m,m=0.8。圖中實(shí)線代表正響應(yīng),虛線代表負(fù)響應(yīng),圓圈代表一維半解析解。

由圖2—圖4可知,本文三維正演計(jì)算的電磁響應(yīng)與一維半解析解吻合很好,僅在符號(hào)反轉(zhuǎn)區(qū)間內(nèi)誤差較大,其余區(qū)間的相對(duì)誤差均在5%以內(nèi)。通過分析發(fā)現(xiàn),在符號(hào)反轉(zhuǎn)區(qū)間內(nèi),電磁響應(yīng)曲線數(shù)值急劇變化,此時(shí)無論是本文算法還是半解析解均無法獲得精確結(jié)果。由此,我們驗(yàn)證了本文的三維正演算法的準(zhǔn)確性。

2.2 典型地質(zhì)體電磁響應(yīng)

對(duì)典型地質(zhì)體模型響應(yīng)進(jìn)行模擬,分別研究半空間、極化地質(zhì)體、非極化地質(zhì)體的響應(yīng)特征,以達(dá)到識(shí)別地下地質(zhì)體激電效應(yīng)特征的目的,為航空電磁數(shù)據(jù)解釋過程中的異常拾取和識(shí)別提供參考。為此,設(shè)計(jì)了如圖5所示的水平板狀體模型。水平板狀體的尺寸為100 m×100 m×30 m,頂部埋深為30 m。飛行系統(tǒng)參數(shù)同上。

首先分別計(jì)算背景電導(dǎo)率為0.010 S/m的均勻非極化半空間、其中嵌入電導(dǎo)率為0.100 S/m的非極化板狀體以及電導(dǎo)率為0.100 S/m的極化板狀體的電磁響應(yīng)(m=0.5,τ=0.01 s,圖6a)。進(jìn)而分別計(jì)算背景電導(dǎo)率為0.010 S/m的均勻極化半空間、其中嵌入電導(dǎo)率為0.100 S/m的非極化板狀體以及電導(dǎo)率為0.100 S/m的極化板狀體的電磁響應(yīng)(m=0.5,τ=0.01 s,圖6b)。

a. 非極化半空間;b. 極化半空間。m=0.5,τ=0.01 s。圖中實(shí)線代表正響應(yīng),虛線代表負(fù)響應(yīng)。

由圖6a可以看出,在t=0.1 ms附近,極化地質(zhì)體響應(yīng)曲線衰減速度比非極化地質(zhì)體慢,兩者衰減速度均較非極化半空間慢,且晚期極化地質(zhì)體曲線出現(xiàn)負(fù)響應(yīng)。研究表明,在含有激電效應(yīng)的地下介質(zhì)中,時(shí)間域電磁響應(yīng)可看成感應(yīng)響應(yīng)和極化響應(yīng)之和[12]。分析可知,在t=0.1 ms附近,低阻地質(zhì)體的存在導(dǎo)致電磁響應(yīng)衰減變慢,由于極化地質(zhì)體產(chǎn)生激發(fā)極化效應(yīng),且極化響應(yīng)與感應(yīng)響應(yīng)同向,因此電磁響應(yīng)衰減速度比非極化地質(zhì)體更慢。由于極化地質(zhì)體的存在,晚期會(huì)出現(xiàn)負(fù)響應(yīng)。

由圖6b可以看出,在t=0.1 ms附近,極化地質(zhì)體響應(yīng)曲線衰減速度比非極化地質(zhì)體慢,兩者衰減速度均較極化半空間慢,在晚期均出現(xiàn)負(fù)響應(yīng),且極化地質(zhì)體負(fù)響應(yīng)出現(xiàn)時(shí)間更早、幅值更大。由于極化圍巖的存在,在晚期3條響應(yīng)曲線均出現(xiàn)符號(hào)反轉(zhuǎn),而極化地質(zhì)體的存在會(huì)使負(fù)值響應(yīng)出現(xiàn)的時(shí)間更早、負(fù)值響應(yīng)幅值更大。

根據(jù)曲線的這些特征差異,我們可以很好地識(shí)別地質(zhì)體是否為極化體,為航空電磁數(shù)據(jù)處理中的異常拾取和識(shí)別提供參考。

2.3 極化參數(shù)對(duì)電磁響應(yīng)的影響特征

為檢驗(yàn)地質(zhì)體電性參數(shù)對(duì)電磁響應(yīng)的影響特征,本文建立一個(gè)三維地質(zhì)體模型,討論其極化參數(shù)對(duì)電磁響應(yīng)的影響。電導(dǎo)率和充電率是影響電磁響應(yīng)的主要參數(shù),本文圍繞這2個(gè)參數(shù)設(shè)計(jì)三維地質(zhì)體模型。地質(zhì)體尺寸為1 000 m×1 000 m×500 m,頂部埋深為100 m(圖7)。

圖7 三維極化地質(zhì)體模型

首先研究極化體電導(dǎo)率對(duì)電磁響應(yīng)的影響特征。假設(shè)極化體電導(dǎo)率分別為1.000、0.200、0.100和0.020 S/m,m=0.5,τ=0.01 s,空氣電導(dǎo)率為10-8S/m,無極化圍巖介質(zhì)電導(dǎo)率為0.010 S/m。圖8展示了不同電導(dǎo)率極化體時(shí)間域航空電磁響應(yīng)。

m=0.5,τ=0.01 s。圖中實(shí)線代表正響應(yīng),虛線代表負(fù)響應(yīng)。

由圖8中給出的對(duì)比結(jié)果可以看出,地質(zhì)體電導(dǎo)率越小,晚期負(fù)響應(yīng)出現(xiàn)越早。這是因?yàn)榈刭|(zhì)體電導(dǎo)率越小,電磁波衰減越快,而由激發(fā)極化效應(yīng)產(chǎn)生的負(fù)響應(yīng)出現(xiàn)越早。

進(jìn)而,我們研究極化體充電率對(duì)電磁響應(yīng)的影響特征。假設(shè)極化體充電率分別為0.0、0.2、0.5、0.8,σ0=0.100 S/m,τ=0.01 s,空氣電導(dǎo)率為10-8S/m,無極化圍巖介質(zhì)電導(dǎo)率為0.010 S/m。不同充電率極化體的電磁響應(yīng)曲線如圖9所示。

σ0=0.100 S/m,τ=0.01 s。圖中實(shí)線代表正響應(yīng),虛線代表負(fù)響應(yīng)。

由圖9可以看出,除了非極化介質(zhì)(m=0.0),其余3種模型均出現(xiàn)了電磁響應(yīng)符號(hào)反轉(zhuǎn)現(xiàn)象。地質(zhì)體充電率越大,負(fù)值響應(yīng)出現(xiàn)越早,且幅值越大。這是因?yàn)闃O化地質(zhì)體充電率越大,激電效應(yīng)越強(qiáng),負(fù)值響應(yīng)出現(xiàn)的時(shí)間也就越早。

3 結(jié)論

1)本文成功提出一種基于非結(jié)構(gòu)有限元的時(shí)間域航空電磁激電效應(yīng)三維正演方法。該方法針對(duì)極化介質(zhì)將Debye模型引入到歐姆定律中,并應(yīng)用二階后推歐拉離散格式對(duì)時(shí)間導(dǎo)數(shù)進(jìn)行離散,實(shí)現(xiàn)電磁場(chǎng)求解。

2)利用不同激發(fā)極化參數(shù)的半空間模型,驗(yàn)證了本文算法具有較高的計(jì)算精度。

3)計(jì)算了半空間和三維地質(zhì)體具有不同的極化條件時(shí)的電磁響應(yīng),發(fā)現(xiàn)可以根據(jù)航空電磁響應(yīng)曲線的衰減特征,有效地識(shí)別圍巖及異常體的極化特征。

4)激電參數(shù)對(duì)時(shí)間域航空電磁系統(tǒng)的影響特征研究表明,地質(zhì)體電阻率越高或者充電率越強(qiáng),極化效應(yīng)影響越嚴(yán)重。

期待本研究可為航空電磁的數(shù)據(jù)處理過程中的異常拾取和識(shí)別提供參考。