層狀多孔介質(zhì)的地震電磁異常數(shù)值模擬方法研究*

于子葉

(中國北京100081中國地震局地球物理研究所)

?

層狀多孔介質(zhì)的地震電磁異常數(shù)值模擬方法研究*

于子葉*

(中國北京100081中國地震局地球物理研究所)

本文針對電激發(fā)模型在各向同性分層介質(zhì)中利用傳播矩陣方法對電磁異常產(chǎn)生過程的數(shù)值模擬． 層狀模型中一部分為飽和多孔介質(zhì)， 另一部分為普通線彈性介質(zhì)， 數(shù)值模擬過程中將固體場作為主要場單獨進行計算， 而將其形變作為等效流體源計算耦合電磁場． 模擬結果顯示， 地震波在經(jīng)過飽和多孔介質(zhì)邊界的時候會產(chǎn)生電磁場， 而其在流體層內(nèi)部傳播的過程中則未產(chǎn)生電磁異常． 模擬結果可以用于解釋一些伴隨地震產(chǎn)生的電磁異常， 包括地震發(fā)生時以及地震波到達時所產(chǎn)生的電磁異?，F(xiàn)象， 另一方面對于前震及其它一些無感地震所產(chǎn)生的電磁異常的解釋也具有一定的參考價值．

電激發(fā)效應 電磁異常 飽和多孔介質(zhì) 電磁場模擬

引言

大量觀測研究表明， 在地震發(fā)生過程中會伴隨產(chǎn)生電磁異常(Garambois， Dietrich， 2001; Honkuraetal， 2002; Matsushimaetal， 2002; 湯吉等， 2008， 2010； Okuboetal， 2011； Frenkel， 2014)． Honkura等(2002)的研究顯示地震波到達之前會產(chǎn)生電磁異常信號， 而Matsushima等(2002)則記錄到電磁異常伴隨地震波產(chǎn)生; Honkura等(2002)則認為地震發(fā)生時檢測到的電磁異常產(chǎn)生于震源， 也就是說， 由于地震破裂而產(chǎn)生了一些電磁場擾動， 而與地震波同時到達的電磁異常則與當?shù)貓龅貤l件有關． 為了解釋這些異常發(fā)展了一些電磁場耦合理論， 例如電激發(fā)理論(Pride， 1994)、 運動電效應理論(Matsushimaetal， 2002)和壓磁效應(Okuboetal， 2011)等． 迄今為止， 地震電磁異常方面的研究多集中于觀測層面， 相關理論及數(shù)值模擬方法仍不甚成熟， 但是， 地震電磁異常對于地震預警乃至地震前兆異常研究起著非常積極的作用(安張輝等， 2013)， 因此一直是地震學研究的熱點問題之一．

介質(zhì)中由于流體作用產(chǎn)生電磁場的機制是頗有前景的研究方向． 巖石實驗研究已表明， 地震波通過含水孔隙介質(zhì)時可產(chǎn)生電磁信號(Zhuetal， 2000； Thompson， Gist， 2012)． 在理論模擬方面， 考慮到層狀介質(zhì)中格林函數(shù)的模擬方法目前已相對成熟， 而且Chen(1999)發(fā)展的反透射系數(shù)法同樣可以高效地計算層狀介質(zhì)中的地震波， 這為層狀介質(zhì)中電磁異常的模擬計算提供了重要的參考． Pride(1994)具體闡述了電激發(fā)原理并給出了約束方程． 基于電激發(fā)原理產(chǎn)生電磁波的過程大致為： 固體場運動使其中的流體產(chǎn)生相對運動(Kennett， 1984)， 而固體顆粒由于物理或化學作用在其上附著電荷， 同時為保持整體的電中性， 流體中攜帶有相反電荷， 流體和固體帶電層的相對運動在宏觀上產(chǎn)生了電磁異常． 隨后， Haartsen和Pride(1997)在約束方程的基礎上給出了全空間電磁格林函數(shù)的表達式； Gao和Hu(2010)以及Ren等(2010)則運用有限差分方法模擬了異常的產(chǎn)生過程； 同時電激發(fā)理論被一些實驗室觀測(Garambois， Dietrich， 2001； Zysermanetal， 2010)所證實． Gao等(2013)根據(jù)傳播矩陣方法改進了積分方式， 提出了支點積分方法． 上述數(shù)值模擬所選取的模型均為孔隙介質(zhì)模型(張丹等， 2013)， 并且震源對于電磁信號的產(chǎn)生具有一定的影響(Hu， Gao， 2011)， 然而由于受數(shù)值模擬方法所限， 對于地震波通過含水層的情況則少有涉及． 鑒于此， 本文擬采用離散波數(shù)的方式對電磁異常信號的產(chǎn)生過程進行數(shù)值模擬， 以期能夠計算不同地層條件下的電磁異常．

本文的模擬實驗基于Pride(1994)的電激發(fā)理論， 不同于前人所使用的將流體場和電磁場作為總體的傳播矩陣進行模擬， 而是在模擬過程中將固體場作為非耦合部分獨立進行模擬計算， 并將已得固體場形變作為流體運動的力源． 震源與孔隙介質(zhì)層分離之后， 震源對于整個電磁異常模擬的作用僅限于產(chǎn)生一個能夠影響到固體場擾動的變量， 這樣在計算過程中即可對一些更復雜的含水地層所產(chǎn)生的電磁異常信號進行數(shù)值模擬．

1 方法原理

1.1 約束方程及固體場計算

圖1b給出了加入含水層的孔隙介質(zhì)模型示意圖， 在該圖中建立柱坐標系， 原點位于地表． 為了使計算簡便， 使z軸通過點源， 向下為正， 第j層邊界記為zj， 地面處記為z0． 根據(jù)Pride(1994)提出的電磁異常約束方程， 在忽略流體場和電磁場的情況下， 介質(zhì)條件簡化為層狀介質(zhì)(圖1a)， 此時約束方程(Haartsen， Pride， 1994， 1997)可簡化為

圖1 普通層狀介質(zhì)模型(a)和加入含水層的孔隙介質(zhì)模型(b)Fig.1 Common stratified media model (a) and stratified media model with fluid-saturated porous layer (b)

(1)

式中，H和G分別為體變模量和剪切模量，ρ為介質(zhì)密度， I為單位矩陣， u為固體場位移， ω為頻率， F為震源， 震源存在于s層中， δjs表示僅在s層中存在震源．

引入矢量柱面諧函數(shù)作為基函數(shù)， 其表達式為

(2)

對固體場向量展開可得(Chen， 1999)

(3)

式中， V為水平和垂直方向上位移和應力的展開系數(shù)，M為系數(shù)矩陣， 依次為

(4)

式(4)的解為

(5)

式中，j為地層編號，Λ為對角矩陣．

結合層間邊界條件和自由邊界條件即可求得每層的固體場方程系數(shù)aj(Chen， 1999).

1.2 耦合電磁場計算

將流體電磁耦合場部分方程(Pride， 1994)整理成兩個關于固體場u， 流體場w和電場E的線性微分方程：

(6)

其中： ρd=iη/(ωκ)， κ為動態(tài)滲透率； εd=ε+i/(ωσ)-ρdL2， C為體應變模量， ρf為液體密度， L為低頻電激發(fā)耦合因子．

由于不存在流體源和電場源， 式(6)中的[C+ωρfI]·u可以作為產(chǎn)生流體運動和電磁場的源． 耦合場的邊界條件為

(7)

接下來， 利用基函數(shù)式(2)分別對電場和流體場進行展開， 可得

(8)

電磁場計算時已經(jīng)將u作為已知源， 即固體場波數(shù)和式(5)的系數(shù)aj是已知的， 將式(8)帶入式(6)整理可得

(9)

取

式(9)是τe的一個關于深度z的四次微分方程， 解為

(10)

其中：

Aj是式(14)所對應的四次齊次微分方程通解的系數(shù)． 求得τe后， 可以推導出用其表示的電場波數(shù)為

(11)

由電場邊界條件及流體場邊界條件式(7)得式(9)的邊界條件為

(12)

利用邊界條件求解可得到τe的4個系數(shù)， 至此流體層中產(chǎn)生電磁場的計算完成．

流體層內(nèi)的系數(shù)可根據(jù)上述方程直接給出， 而流體層外的數(shù)據(jù)則需根據(jù)電磁場邊界條件繼續(xù)計算. 在不存在流體作用和自由電荷影響的情況下， 假設其它層為各項同性均勻介質(zhì)， 電場方程則可簡化為

(13)

其中

式中σ為電導率，μ為磁導率，ε為介質(zhì)的介電常數(shù)， 其邊界條件為

(14)

利用式(8)中電場等式展開后， 帶入麥克斯韋方程可得無流體時電場系數(shù)方程為

(15)

求解得

(16)

式中： C為方程系數(shù)， 通過邊界條件即可確定； η為式(15)所對應二次方程的根．

至此， 固體場、 流體場、 電磁場波數(shù)域展開系數(shù)已全部求出， 對其反變換后可得其頻率域的空間解為

(17)

圖2 電磁異常模擬程序結構框圖Fig.2 Layout of EM anomaly simulation program

式中，θ為觀測點坐標與x軸的夾角． 固體場與電場矢量反變換形式相同， 均用v表示， 對v進行傅里葉反變換得到其時間域的表達式為

(18)

至此頻率域矢量求解完畢．

1.3 程序結構

本文的電磁異常數(shù)值計算過程如圖2所示． 首先， 利用類傳播矩陣方法計算固體場波數(shù)域系數(shù)， 結合所讀取的含水層參數(shù)計算電磁場波數(shù)域系數(shù)； 然后， 根據(jù)所得到的電磁場進行反變換求得頻率域解； 最后， 進行傅里葉反變換求得其時間域解．

由于Λj存在一個極小的數(shù)， 若在系數(shù)a的求解過程中將其當成一個整體矩陣求解， 則很容易出現(xiàn)病態(tài)矩陣， 故本文采用Chen(1999)的計算方式， 即將整體矩陣分解成小的傳播矩陣， 從而降低矩陣的病態(tài)性， 但這樣使得式(5)中系數(shù)a的計算過程相對復雜了一些．

模擬程序采用標準C(C99標準)進行編寫， 利用Window平臺下的MinGW編譯器(GCC版本4.7.0)進行編譯， 便于該程序移植到Linux平臺．

2 模擬計算與結果比較

2.1 震源模型

本文在研究過程中著重于固體場擾動產(chǎn)生電磁波的過程， 所以震源的作用在于產(chǎn)生地震波， 從而能影響到相關的流體場， 而震源本身的形式并不重要， 因此在模擬中采用相對簡單的矩張量點源模型， 震源深度為15 km， 其表達式為

(19)

其中，

(20)

式中：M0=1.5×108N·m， 相當于M6.0地震； 震源時間函數(shù)s(t)為地震子波， 表達式為

(21)

式中，Tc=4 s，f0=1 Hz.

2.2 地層模型

模擬的地層為各向同性介質(zhì)， 共4層． 本文著重研究地震波傳播經(jīng)過層狀介質(zhì)時所產(chǎn)生的電磁波， 為了避免地層間的多次反射波干擾， 將地層參數(shù)設為同一數(shù)值． 所選取的彈性參數(shù)均為典型巖石的彈性參數(shù)． 其中2號地層為含水層， 其參數(shù)列于表1． 表中孔隙度和鹽濃度的選取參考了Gao等(2013)的模擬參數(shù)， 電導率σ0和低頻電激發(fā)耦合因子L0為鹽濃度C0與電勢ζ的函數(shù)， 其中電勢ζ為飽和的石英介質(zhì)25℃時在pH=7的條件下所得的實驗值(Haartsen， Pride， 1997)．

表1 含水層參數(shù)

2.3 數(shù)值結果

將接收點設在含水層界面上， 該點距震中10 km， 方位角為30°， 測定地震波振幅和電場強度． 觀測點之所以這樣選取， 著重于觀測地震波在入射和出射時所產(chǎn)生的電磁異常， 計算結果如圖3和4所示．

數(shù)值模擬所得圖形采用歸一化表示， 再乘以相應的系數(shù)即所得電場大?。?從圖3與圖4的對比可以看出: 初始接收到的電場與固體波場的到時相差約1.7 s， 這是由于地震波與電磁波的速度存在差異所致; 地震波進入含水層下界面時產(chǎn)生了電磁場， 二者從下層界面同時出發(fā)， 但電磁場提前到達； 電磁場信號與地震波信號呈現(xiàn)出很強的相關性， 而且僅在地震波入射和出射含水層時產(chǎn)生電磁異常； 地震波振幅與所產(chǎn)生的相應的電場振幅呈正相關， 且電磁波較地震波衰減得快．

圖3 含水層上界面固體場位移

圖4 含水層上界面電場

2.4 觀測結果比較

現(xiàn)有很多研究均報道直接觀測到了同震電磁信號， 其中比較有代表性的是2013年甘肅岷縣漳縣MS6.6地震的同震電磁信號(楊皓等， 2015)． 該地震中， 電場異常與地震波同時到達， 且形狀相似． 當P波到達時， 電場的水平和垂直方向均產(chǎn)生了明顯的電磁異常， 而之后S波到達時所產(chǎn)生的電磁異常幅度較P波所產(chǎn)生的要大， 這與地震波振幅相關． 假設接收點以下介質(zhì)均為飽和含水介質(zhì)， 那么觀測結果與數(shù)值模擬得到的結果較為一致， 即在地震波入射和出射含水層界面時產(chǎn)生電磁異常， 且其異常幅度與地震波振幅相關； 但電場曲線與地震動曲線相比缺少高頻特征， 這說明在電磁場與波場的耦合過程中高頻部分耦合得不好， 可能與流體性質(zhì)有關， 在高頻運動過程中由于流動性等因素導致流體場與固體場的相對運動減弱， 使得所產(chǎn)生的電磁場減弱， 這也說明Pride(1994)的約束方程可能缺少一個高頻阻尼項． 另外， 未觀測到與地震破裂過程伴隨產(chǎn)生的電磁異常， 這可能是電磁場隨深度衰減過快所致， 也可能是在震源破裂過程中并未產(chǎn)生電磁波．

3 討論與結論

由式(11)可知， 電場波數(shù)完全取決于固體場波數(shù)域展開系數(shù)， 而固體場的波數(shù)解則取決于震源， 即電磁場的解是震源的函數(shù)， 所以震源函數(shù)與電磁場信號存在一定的相關性， 這說明地電磁異常中常用電磁信號與地震信號進行互相關計算是合理的． 式(11)轉換成通用形式如下：

(22)

式中，M和N為系數(shù)． 可以看出， 電磁場與地震波場存在著一定的相關性． 整理后，Mi中包含了ω0,ω1,ω2和ω3項. 由傅里葉變換的導數(shù)性質(zhì)可知， 電磁場還與固體場的時間域?qū)?shù)相關．

由模擬結果可以看到兩個獨立的信號， 這說明電磁場僅在震動通過電性變化區(qū)域時才會產(chǎn)生， 這也說明如果在實際中由電激發(fā)效應產(chǎn)生電磁場， 必要條件是介質(zhì)孔隙層間流體條件不一致． 在前震、 地震成核階段以及無感地震中， 均會有一定程度的地下巖石局部形變， 若該形變場通過孔隙條件復雜的流體層， 那么由數(shù)值模擬結果可知該過程會產(chǎn)生電磁場， 這就為解釋一些同震電磁異常以及無感地震中觀測到的電磁場異常提供了一種可能的解釋． 鑒于地下條件的復雜性， 電磁異常的產(chǎn)生過程可能并非如Pride(1994)電激發(fā)理論中所描述的那樣簡單， 本研究對地震波場、 流體場和電磁場等3場耦合產(chǎn)生電磁場的過程進行了有現(xiàn)實意義的探索．

本文僅計算了最強的一次場， 并未考慮電磁場在地層中的反射以及不同頻率電磁波的衰減， 希望在以后研究中予以改進． 另外， 本文僅用簡化的點源進行了解釋和模擬， 而未用實際震源進行更具體的模擬， 進一步研究中需要帶入真實的震源模型， 并與實際地震進行對比．

安張輝， 杜學彬， 譚大誠， 范瑩瑩， 劉君， 崔騰發(fā). 2013. 四川蘆山MS7.0和汶川MS8.0地震前地電場變化研究[J]. 地球物理學報， 56(11): 3868--3876.

An Z H， Du X B， Tan D C， Fan Y Y， Liu J， Cui T F. 2013. Study on the geo-electric field variation of Sichuan LushanMS7.0 and WenchuanMS8.0 earthquake[J].ChineseJournalofGeophysics， 56(11): 3868--3876 (in Chinese).

湯吉， 詹艷， 王立鳳， 徐建郎， 趙國澤， 陳小斌， 董澤義， 肖騎彬， 王繼軍， 蔡軍濤， 徐光晶. 2008. 5月12日汶川8.0級地震強余震觀測的電磁同震效應[J]. 地震地質(zhì)， 30(3): 739--745.

Tang J， Zhan Y， Wang L F， Xu J L， Zhao G Z， Chen X B， Dong Z Y， Xiao Q B， Wang J J， Cai J T， Xu G J. 2008. Coseismic signal associated with aftershock of theMS8.0 Wenchuan earthquake[J].SeismologyandGeology， 30(3): 739--745 (in Chinese).

湯吉， 詹艷， 王立鳳， 董澤義， 趙國澤， 徐建郎. 2010. 汶川地震強余震的電磁同震效應[J]. 地球物理學報， 53(3): 526--534.

Tang J， Zhan Y， Wang L F， Dong Z Y， Zhao G Z， Xu J L. 2010. Electromagnetic coseismic effect associated with aftershock of WenchuanMS8.0 earthquake[J].ChineseJournalofGeophysics， 53(3): 526--534 (in Chinese).

楊皓， 詹艷， 趙凌強， 陳小斌， 姜峰. 2015. 2013年甘肅岷縣漳縣MS6.6地震震電磁信號和同震信號觀測[J]. 震災防御技術， 10(增刊): 785--793.

Yang H， Zhan Y， Zhao L Q， Chen X B， Jiang F. 2015. The seismo-electromagnetic signal and coseismic signal of 2013 Minxian-ZhangxianMS6.6 earthquake[J].TechnologyforEarthquakeDisasterPrevention， 10(Suppl): 785--793 (in Chinese).

張丹， 任恒鑫， 黃清華. 2013. 孔隙介質(zhì)地震電磁信號的數(shù)值模擬研究[J]. 地球物理學報， 56(8): 2739--2747.

Zhang D， Ren H X， Huang Q H. 2013. Numerical simulation study of co-seismic electromagnetic signals in porous media[J].ChineseJournalofGeophysics， 56(8): 2739--2747 (in Chinese).

Chen X F. 1999. Seismogram synthesis in multi-layered half-space: PartⅠ.Theoretical formulations[J].EarthquakeResearchinChina， 13(2): 149--174.

Frenkel J. 2014. On the theory of seismic and seismoelectric phenomena in a moist soil[J].JEngMech， 131(9): 879--887.

Gao Y X， Hu H S. 2010. Seismoelectromagnetic waves radiated by a double couple source in a saturated porous medium[J].GeophysJInt， 181(2): 873--896.

Gao Y X， Chen X F， Hu H S， Zhang J. 2013. Early electromagnetic waves from earthquake rupturing:Ⅰ.Theoretical formulations[J].GeophysJInt， 192(3): 1288--1307.

Garambois S， Dietrich M. 2001. Seismoelectric wave conversions in porous media: Field measurements and transfer function analysis[J].Geophysics， 66(5): 1417--1430.

Haartsen M W， Pride R S. 1994. Modeling of coupled electroseismic wave propagation from point sources in layered media[C]∥SEGTechnicalProgramExpandedAbstracts. Los Angeles, California: SEG, 13(1): 1155--1158.

Haartsen M W， Pride S R. 1997. Electroseismic waves from point sources in layered media[J].JGeophysRes， 102(B11): 24745--24769.

Honkura Y， Matshshima M， Oshiman N， Tun?er M K， Bari， Ito A， Iio Y， Iikara A M. 2002. Small electric and magnetic signals observed before the arrival of seismic wave[J].EarthPlanetsSpace， 54(12): e9--e12.

Hu H S， Gao Y X. 2011. Electromagnetic field generated by a finite fault due to electrokinetic effect[J].JGeophysRes， 116(B8): B08302.

Kennett B L N. 1984. Seismic wave propagation in stratified media[J].GeophysJRastrSoc， 86(1): 219--220.

Matsushima M， Honkura Y， Oshiman N， Baris S， Tuncer M K， Tank S B， Celik C， Takahashi F， Nakanishi M， Yoshimura R， Petras R， Komut T， Tolak E， Ito A， Iio Y， Isikara A M. 2002. Seismo-electromagnetic effect associated with the Izmit earthquake and its aftershocks[J].BullSeismolSocAm， 92(1): 350--360.

Okubo K， Takeuchi N， Utsugi M， Yumoto K， Sasai Y. 2011. Direct magnetic signals from earthquake rupturing: Iwate-Miyagi earthquake ofM7.2， Japan[J].EarthPlanetSciLett， 305(1/2): 65--72.

Pride S R. 1994. Governing equations for the coupled electromagnetics and acoustics of porous media[J].PhysRevB， 50(21): 15678--15696.

Ren H X， Huang Q H， Chen X F. 2010. A new numerical technique for simulating the coupled seismic and electromagnetic waves in layered porous media[J].EarthquakeScience， 23(2): 167--176.

Thompson A H， Gist G A. 2012. Geophysical applications of electrokinetic conversion[J].LeadEdge， 12(12): 1169--1173.

Zhu Z Y， Haartsen M W， Toks?z M N. 2000. Experimental studies of seismoelectric conversions in fluid-saturated porous media[J].JGeophysRes， 105(B12): 28055--28064.

Zyserman F I， Gauzellino P M， Santos J E. 2010. Finite element modeling of SHTE and PSVTM electroseismics[J].JApplGeophys， 72(2): 79--91.

On numerical simulation method of the seismic electromagnetic anomaly in porous stratified media

Yu Ziye*

(InstituteofGeophysics，ChinaEarthquakeAdministration，Beijing100081，China)

This paper focuses on numerical simulation of the seismic electro-magnetic anomaly generation process in isotropic layered medium by means of propagation-matrix method. For the layered model， some are fluid-saturated porous medium， and the others are linear elastic medium. During the simulation， the solid field is taken as a principal field and calculated individually， and its deformation is taken as the equivalent fluid source. The simulation results show that when the seismic wave goes through the boundary of the fluid-saturated porous medium， the electromagnetic (EM) wave is generated simultaneously， which can be used to illustrate the seismic EM anomaly accompanied by seismic source or the seismic wave. On the other hand， it can elucidate EM signal coupled with the foreshock or the tremor.

electrokinetic effect； electromagnetic anomaly； fluid-saturated porous medium； electromagnetic field simulation

國家重大科學基礎設施項目東半球空間環(huán)境地基綜合監(jiān)測子午鏈“子午工程”資助.

2016-01-06收到初稿， 2016-04-13決定采用修改稿.

10.11939/jass.2016.06.006

P319.1+2

A

于子葉. 2016. 層狀多孔介質(zhì)的地震電磁異常數(shù)值模擬方法研究. 地震學報, 38(6): 869--877. doi:10.11939/jass.2016.06.006.

Yu Z Y. 2016. On numerical simulation method of the seismic electromagnetic anomaly in porous stratified media.ActaSeismologicaSinica, 38(6): 869--877. doi:10.11939/jass.2016.06.006.

*通訊作者 e-mail: cangye@hotmail.com