基于雙電層模型的頻域震電響應(yīng)數(shù)值模擬

韓 麗, 胥鐵瀟, 金勝昔, 董文宇, 嵇艷鞠

吉林大學(xué)儀器科學(xué)與電氣工程學(xué)院，長春 130026

0 引言

地震的孕育和發(fā)生過程中常常會出現(xiàn)電磁異常擾動。研究認(rèn)為地震波到達(dá)之前產(chǎn)生的電磁信號是與震源相關(guān)的，而與地震波同時到達(dá)的電磁信號則與當(dāng)?shù)氐膱龅貤l件相關(guān)[1-2]。但目前對這些地震電磁信號產(chǎn)生的機(jī)制仍存有爭議。

許多學(xué)者推測雙電層動電效應(yīng)是最有可能的機(jī)制之一，這是因?yàn)閯与娦?yīng)要求介質(zhì)是孔隙介質(zhì)，是普遍存在于地殼表層地層中的[3-7]。動電效應(yīng)與孔隙流體的運(yùn)動是密切相關(guān)的，通過地震孕育和發(fā)生過程中的動電效應(yīng)研究，進(jìn)一步了解地震中引發(fā)電磁異?，F(xiàn)象的機(jī)制，對于地震預(yù)報和監(jiān)測、防震減災(zāi)都有重要意義[8-9]。

傳統(tǒng)的震電響應(yīng)建模中，電磁部分以確定的地震場為源，采用準(zhǔn)靜態(tài)近似的方式，通過求解電勢的泊松方程獲得電場[10-17]。這種準(zhǔn)靜態(tài)近似的方式可以降低模擬難度，能一定程度提高計(jì)算效率；然而準(zhǔn)靜態(tài)近似在震電響應(yīng)模擬中會帶來數(shù)值誤差，尤其是在高礦化度條件下，誤差特別顯著[4]。并且這種方法必須滿足準(zhǔn)靜態(tài)近似條件，即模型區(qū)域必須在靜電近場范圍內(nèi)，頻率也位于一個較低的區(qū)間范圍。

本文基于Maxwell全方程的頻域有限差分方法，進(jìn)行了震電響應(yīng)的頻率電磁響應(yīng)數(shù)值模擬，并分析和討論在均勻彈性介質(zhì)中，電導(dǎo)率層狀和電導(dǎo)率異常體模型下的震電特征，其中模型均沒考慮自由地表影響。

1 基于Maxwell全方程的頻域震電響應(yīng)模擬

有限差分方法是采用數(shù)值近似的方法，將原本連續(xù)的函數(shù)離散化，形成線性方程組，之后借助計(jì)算機(jī)來進(jìn)行大規(guī)模計(jì)算的一種方法。

本次借助頻域有限差分方法，可以大大簡化在復(fù)雜空間中進(jìn)行波場計(jì)算的計(jì)算量，在二維空間中將空間分割成一個個網(wǎng)格(圖1)。

Ex1--Ex6為網(wǎng)格棱邊上x方向的電場分量；Hz1--Hz4為網(wǎng)格面上z方向的磁場分量。

圖1 二維頻域有限差分圖

Fig1 Two dimensional FDFD

在計(jì)算時每一個網(wǎng)格單獨(dú)計(jì)算，這樣在計(jì)算復(fù)雜空間模型時，不會因復(fù)雜的邊界條件而大大增加計(jì)算成本。地震波造成動電效應(yīng)產(chǎn)生電磁波頻率較低，此時可以忽略位移電流，將Maxwell方程簡化成：

×E=-iωμ0H;

(1)

×H=σE;

(2)

·E=0;

(3)

·H=0。

(4)

式中：ω為角頻率；E為垂直極化的橫波(PSV)和水平極化的磁場(TM)波產(chǎn)生的電場；H為磁場；i為虛數(shù)單位。其他參數(shù)參見表1。

與單一電磁場頻域有限差分不同的是，在動電效應(yīng)中，產(chǎn)生的電磁波是由地震波在介質(zhì)中傳播所激發(fā)的，因此在場中沒有電流源，滲流位移代替電流成為了電磁場的初始條件之一。

本次模擬耦合地震波和電磁(EM)波的理論基礎(chǔ)是Pride震電耦合控制方程。為簡化方程，我們忽略了地震感應(yīng)電磁場的反饋，在前半部分以格林函數(shù)解析的方式來求解地震B(yǎng)iot方程；而對于電磁部分，不同于傳統(tǒng)引入準(zhǔn)靜態(tài)近似的方式，我們將采用Maxwell全方程來求解電磁波場。

表1 介質(zhì)參數(shù)、符號、值和SI單位

考慮二維有限域各向同性的流體飽和多孔介質(zhì)，通過丟棄電過濾反饋，Pride震電耦合控制如下：

(5)

(6)

(7)

(8)

(9)

×E=-μH。

(10)

式中：v為固相粒子速度；q為平均相對流體速度；τ為應(yīng)力；ρ=(1-φ)ρs+φρf，為體積密度；ρm=Tρf/φ,為有效的流體密度；P為孔隙流體壓力；i和j為變量x或者y的方向；以固相粒子速度v為例，vi,j為vi在j方向的空間偏導(dǎo)數(shù)，vj,ji為變量vj在j方向和i方向的空間偏導(dǎo)數(shù)之和；變量上方·代表該變量對時間的導(dǎo)數(shù)；α、λ和M為彈性模量。α、λ和M由下式給出：

α=1-Kb/Ks；

(11)

(12)

(13)

其余波場參數(shù)意義及設(shè)置如表1所示。

首先使用格林函數(shù)方式求解獨(dú)立的地震波場結(jié)果[16]。其中心思想是在三維空間傅里葉時間-拉普拉斯變換域中建立方程，場量的F(x,t)的時間-拉普拉斯變換：

(14)

式中：t為時間；s為拉普拉斯變換參數(shù)，可以選擇s=iω將其簡化為時間傅里葉變換形式。則其三維空間傅里葉變換為

(15)

式中：kn為波數(shù)向量k的第n個笛卡爾向量分量；xn為向量x第n個分量的位置坐標(biāo)。

以往的震電模擬中使用時域有限差分算法求解。時域有限差分相比頻域有限差分算法占用內(nèi)存更少，可以計(jì)算大范圍區(qū)域的電磁場，但是每一次迭代都要對所有網(wǎng)格全部刷新一遍，計(jì)算速度較慢。當(dāng)計(jì)算地震波和電磁波耦合的波場時，電磁波和地震波速度上的巨大差異給計(jì)算帶來了很大困難。在空間域上，為避免差分算法誤差的影響，空間波長要依據(jù)地震波波長來確定；在時域上，時間步長又需依據(jù)電磁波速度來確定，兩者巨大的差距會給結(jié)果帶來很大的誤差。因此本文選擇頻域有限差分方法來進(jìn)行計(jì)算。

由解耦狀態(tài)下的電磁場得到電磁控制方程如下：

(16)

(17)

只計(jì)算xOy平面上的波場情況，可以將對z方向上的求導(dǎo)結(jié)果設(shè)置為0，將旋度在直角坐標(biāo)下分解并簡化可得：

(18)

(19)

(20)

式中：qx、qy和Ex、Ey分別為x、y方向的電場分量；Hz、Hy分別為z、y方向的磁場分量。之后再對其進(jìn)行離散，采用中心差分項(xiàng)代替微分項(xiàng)，可以得到：

(21)

(22)

(23)

式中，Δx、Δy分別表示在x、y方向的網(wǎng)格步長。

最后對頻域差分方程求解，可以解出各個頻點(diǎn)的電場值。

2 數(shù)值模擬

模型大小均設(shè)置為700 m×700 m，網(wǎng)格128×128。此外，源定義在中心網(wǎng)格節(jié)點(diǎn)處，接收點(diǎn)定義在網(wǎng)格節(jié)點(diǎn)(90, 90)。在水平應(yīng)力項(xiàng)上加入中心頻率為f0= 25 Hz的瑞雷子波作為震源函數(shù)，根據(jù)格林函數(shù)方式[15]計(jì)算所得地震滲流速度結(jié)果。具體參數(shù)設(shè)置如表1所示。

在如上設(shè)置的全空間模型下，震電響應(yīng)產(chǎn)生的電場模擬結(jié)果如圖2所示，可以明顯觀察出其在全空間模型中的頻域波場解。

為了分析復(fù)雜介質(zhì)條件下的震電耦合特性，我們設(shè)置一個雙層層狀介質(zhì)模型。設(shè)置上半層-180～350 m電導(dǎo)率為1×10-2S/m，下半層-350～-180 m電導(dǎo)率為1×10-3S/m(圖3)，其他介質(zhì)參數(shù)見表1。由圖3可以看出，層狀模型和全空間模型中的電場分布不同，可以明顯看到電場在分界面的反射和折射現(xiàn)象。

為了模擬現(xiàn)實(shí)地質(zhì)條件的復(fù)雜情況，在全空間模型中加入了電導(dǎo)率不同的異常體，大地電導(dǎo)率為1×10-3S/m，異常體電導(dǎo)率為 2 S/m，異常體為設(shè)置在x方向140 ～200 m、y方向140 ～200 m處的正方形(圖4)，其他介質(zhì)參數(shù)見表1。由圖4可以看出，電場在異常體處發(fā)生了明顯的反常擴(kuò)散現(xiàn)象，通過圖像可以清晰地定位異常體的位置。

圖2 全空間模型的電場頻域波場圖

圖3 層狀模型的電場頻域波場圖

圖4 異常體模型的電場頻域波場圖

通過層狀和異常體的模型與全空間模型的對比，可以看出這種頻域有限差分方法可以有效模擬震電響應(yīng)，并且識別不同介質(zhì)分界面。

3 結(jié)語

本文提出了一種基于Maxwell全方程的頻域有限差分方法，模擬分析了以雙電層理論為基礎(chǔ)的震電耦合效應(yīng)。在彈性均勻介質(zhì)中，以電導(dǎo)率雙層層狀和電導(dǎo)率異常體模型為實(shí)驗(yàn)案例，對比全空間模型，分析了震電耦合特性。結(jié)果表明本方法能夠有效地捕獲對比介質(zhì)界面上地震與電磁波場的反射和透射現(xiàn)象，識別異常位置，突出了地震響應(yīng)模擬檢測地下介質(zhì)特性的能力。