基于Debye模型的時間域航空電磁激電效應正演模擬

張鑫崇,殷長春

吉林大學地球探測科學與技術學院,長春 130026

0 引言

時間域航空電磁法采用機載平臺搭載電磁儀器,可以在沙漠、高山、沼澤、湖泊等地形條件復雜的地區進行勘探,是一種高效的地球物理手段,現已被廣泛應用于礦產資源、地下水及環境工程等眾多勘探領域[1]。與頻率域電磁法相比,時間域航空電磁法對深部目標具有更好的分辨率和探測靈敏度[2]。

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

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

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

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

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

1 方法原理

1.1 控制方程

時間域麥克斯韋方程可以寫成:

(1)

(2)

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

(3)

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

(4)

(5)

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

(6)

1.2 極化介質歐姆定律

傳導電流密度可由歐姆定律描述,即

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

(7)

式中,σ為電導率,在非極化介質中是不變的。然而,當存在激發極化效應時,電導率是頻散的,可用Debye模型[24]來描述:

(8)

式中:σ0為零頻電導率;m為充電率;τ為時間常數;i為虛數單位;ω為角頻率。此時,電流密度和電場強度不再滿足式(7)簡單的線性關系。下面推導極化介質中歐姆定律的形式。

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

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

(9)

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

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

(10)

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

(11)

式(11)即為極化介質中的歐姆定律,其左右兩端均出現了時間導數。依然采用二階后推歐拉離散格式對其進行離散,則式(11)可整理為

(12)

為方便后續推導,將式(12)相關系數用一些字符替代并整理為

(13)

其中:

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

(14)

1.3 有限元方法

本文使用矢量有限元法對式(14)進行空間離散。對于非結構四面體網格,第n個網格r處tk時刻的電場強度可表示為

(15)

(16)

式中,j=1, 2, ???, 6。令加權殘差Rn=0,并通過一系列代數運算,可以建立第n個網格的如下有限元方程:

(17)

式中,相關的矩陣可寫為:

(18)

將全部網格的式(18)組裝起來,可以得到tk時刻電場強度的稀疏對稱方程組,即

[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)

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

Ae(tk)=b(tk)。

(20)

式中:A為系數矩陣;e(tk)為模型剖分網格單元各棱邊上tk時刻的待求電場強度;b(tk)為tk時刻的右端源項。

我們使用直接求解器MUMPS求解得到各棱邊的電場[28],而電場和磁場對時間的導數可以分別用矢量基函數及其導數通過插值計算得到。

2 數值算例

2.1 精度驗證

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

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

空氣電導率設為10-8S/m,分別計算不同均勻半空間電導率、充電率、時間常數的電磁響應,并與一維半解析解進行對比(圖2—圖4)。

a. 電磁響應;b. 與一維半解析解的相對誤差。σ0=0.010 S/m,m=0.8。圖中實線代表正響應,虛線代表負響應,圓圈代表一維半解析解。

由圖2—圖4可知,本文三維正演計算的電磁響應與一維半解析解吻合很好,僅在符號反轉區間內誤差較大,其余區間的相對誤差均在5%以內。通過分析發現,在符號反轉區間內,電磁響應曲線數值急劇變化,此時無論是本文算法還是半解析解均無法獲得精確結果。由此,我們驗證了本文的三維正演算法的準確性。

2.2 典型地質體電磁響應

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

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

a. 非極化半空間;b. 極化半空間。m=0.5,τ=0.01 s。圖中實線代表正響應,虛線代表負響應。

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

由圖6b可以看出,在t=0.1 ms附近,極化地質體響應曲線衰減速度比非極化地質體慢,兩者衰減速度均較極化半空間慢,在晚期均出現負響應,且極化地質體負響應出現時間更早、幅值更大。由于極化圍巖的存在,在晚期3條響應曲線均出現符號反轉,而極化地質體的存在會使負值響應出現的時間更早、負值響應幅值更大。

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

2.3 極化參數對電磁響應的影響特征

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

圖7 三維極化地質體模型

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

m=0.5,τ=0.01 s。圖中實線代表正響應,虛線代表負響應。

由圖8中給出的對比結果可以看出,地質體電導率越小,晚期負響應出現越早。這是因為地質體電導率越小,電磁波衰減越快,而由激發極化效應產生的負響應出現越早。

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

σ0=0.100 S/m,τ=0.01 s。圖中實線代表正響應,虛線代表負響應。

由圖9可以看出,除了非極化介質(m=0.0),其余3種模型均出現了電磁響應符號反轉現象。地質體充電率越大,負值響應出現越早,且幅值越大。這是因為極化地質體充電率越大,激電效應越強,負值響應出現的時間也就越早。

3 結論

1)本文成功提出一種基于非結構有限元的時間域航空電磁激電效應三維正演方法。該方法針對極化介質將Debye模型引入到歐姆定律中,并應用二階后推歐拉離散格式對時間導數進行離散,實現電磁場求解。

2)利用不同激發極化參數的半空間模型,驗證了本文算法具有較高的計算精度。

3)計算了半空間和三維地質體具有不同的極化條件時的電磁響應,發現可以根據航空電磁響應曲線的衰減特征,有效地識別圍巖及異常體的極化特征。

4)激電參數對時間域航空電磁系統的影響特征研究表明,地質體電阻率越高或者充電率越強,極化效應影響越嚴重。

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