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

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

0 引言

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

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

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

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

1 基于Maxwell全方程的頻域震電響應模擬

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

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

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

圖1 二維頻域有限差分圖

Fig1 Two dimensional FDFD

在計算時每一個網(wǎ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。

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

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

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

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

(5)

(6)

(7)

(8)

(9)

×E=-μH。

(10)

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

α=1-Kb/Ks；

(11)

(12)

(13)

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

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

(14)

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

(15)

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

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

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

(16)

(17)

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

(18)

(19)

(20)

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

(21)

(22)

(23)

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

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

2 數(shù)值模擬

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

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

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

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

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

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

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

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

3 結(jié)語

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