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

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

0 引言

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

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

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

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

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

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

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

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

圖1 二維頻域有限差分圖

Fig1 Two dimensional FDFD

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

×E=-iωμ0H;

(1)

×H=σE;

(2)

·E=0;

(3)

·H=0。

(4)

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

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

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

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

考慮二維有限域各向同性的流體飽和多孔介質，通過丟棄電過濾反饋，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ù)意義及設置如表1所示。

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

(14)

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

(15)

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

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

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

(16)

(17)

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

(18)

(19)

(20)

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

(21)

(22)

(23)

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

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

2 數(shù)值模擬

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

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

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

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

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

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

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

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

3 結語

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