孔隙介質(zhì)地震電磁信號的數(shù)值模擬研究

張 丹，任恒鑫，黃清華＊

1 北京大學(xué)地球與空間科學(xué)學(xué)院地球物理學(xué)系，北京 100871

2 中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院，合肥 230026

1 引 言

大量野外觀測研究表明，許多地震過程中伴隨發(fā)生不同程度的電磁輻射異?，F(xiàn)象，例如，Matsushim等［1］記錄了1999年I·zmit7.4級地震及其余震過程中的電場和磁場，Karakelian等［2］記錄到了1999年Hector Mine 7.1 級地震過程中的電磁場，湯吉等［3－4］在武都漢王地震臺附近監(jiān)測到了汶川地震余震事件的同震電磁信號，并指出電磁場和地震波在波形和頻率上都具有良好的一致性.事實上早在1939年Ivanov［5］就發(fā)現(xiàn)了流動電勢效應(yīng)——E效應(yīng)，即在不施加外加電場的情況下，記錄到了地震波產(chǎn)生的電場.為解釋這一現(xiàn)象，F(xiàn)renkel［6］建立了一套震電方程組，該方程組描述了孔隙介質(zhì)中固－液相對運動而產(chǎn)生的流動電流，但是并未考慮輻射電磁波的問題.在此基礎(chǔ)上，Pride［7］通過體積平均方法推導(dǎo)出了一套更為完整的流體飽和孔隙介質(zhì)中的震電波場控制方程組，為之后的震電效應(yīng)研究提供了基礎(chǔ).

在震電問題的數(shù)值模擬方面，Haartsen ＆Pride［8］利用全局矩陣方法對層狀介質(zhì)中點源激發(fā)的震電波場進(jìn)行了研究，結(jié)果表明存在兩種類型的電磁波：一種是地震波在其經(jīng)過的地方激發(fā)的伴隨電磁波，其存在依賴于地震波的傳播；另一種是當(dāng)?shù)卣鸩ń?jīng)過介質(zhì)分界面時，會產(chǎn)生獨立于地震波傳播的輻射電磁場，以電磁波速度傳播.Garambois ＆Dietrich［9］進(jìn)一步對水平層狀介質(zhì)中點震源激發(fā)的震電波場進(jìn)行了系統(tǒng)的研究，并探討了電磁信號對介質(zhì)參數(shù)變化的響應(yīng).但整體而言，針對孔隙介質(zhì)震電波場的數(shù)值模擬研究并不太多.最近，隨著合成地震圖的數(shù)值算法的發(fā)展（例如，Chen［10］和 Ge ＆Chen［11］發(fā)展的廣義反透射系數(shù)方法已在合成地震圖的數(shù)值模擬研究中被證實具有穩(wěn)定、高效和適用范圍較廣等特點），Ren等［12－13］將廣義反透射系數(shù)方法拓展到了孔隙介質(zhì)中震電波場的數(shù)值模擬研究中，得到了一套計算點源激發(fā)的震電波場的數(shù)值模擬方法.

在震電的野外試驗和室內(nèi)實驗方面，Martner ＆Sparks［14］，Thompson ＆ Gist［15－16］和 Zhu等［17－19］等工作均表明含流體孔隙介質(zhì)中的地震波可以激發(fā)電磁信號，利用震電轉(zhuǎn)換可以探測一定深度的油氣儲層分界面.由于震電效應(yīng)能同時反映介質(zhì)彈性參數(shù)和電性參數(shù)的變化，基于震電效應(yīng)的震電勘探作為一種潛在的油氣勘探方法受到了地球物理學(xué)界的重視.

上述學(xué)者的研究大多側(cè)重勘探領(lǐng)域（例如，數(shù)值模擬局限于爆炸點源或單力點源），雖然震電效應(yīng)在天然地震監(jiān)測和預(yù)警研究方面也具有潛在的應(yīng)用價值，但是與之相關(guān)的報道和研究甚少［20］.目前，Gao等［21］模擬了雙力偶點源激發(fā)的地震電磁波場，Hu等［22］研究分析了有限斷層在地表附近激發(fā)的地震電磁波場，Ren等［23］將爆炸源和單力點源逐步拓展到雙力偶點源和有限斷層面源的模型，實現(xiàn)了相應(yīng)模型的震電波場數(shù)值模擬算法，Gao等［24－25］模擬了地震發(fā)生時刻產(chǎn)生的電磁信號，并研究了該信號隨介質(zhì)參數(shù)的變化響應(yīng).但以上工作尚未深入研究震源類型等參數(shù)對結(jié)果的影響.對孔隙介質(zhì)中震電波場認(rèn)識比較有限的現(xiàn)狀直接制約了地震電磁機(jī)理研究的進(jìn)一步開展.有鑒于此，本文利用Ren等［23］的方法對雙力偶點源以及有限斷層面源在不同模型下產(chǎn)生的電磁信號進(jìn)行了數(shù)值模擬，考察了地震波和電磁信號的相關(guān)性，并討論了不同的震源時間函數(shù)類型和不同的Ricker子波峰值頻率對地震和電磁信號的影響.

2 方 法

Pride［7］用體積平均方法建立了一套震電波場宏觀控制方程組，該方程組描述了流體飽和孔隙介質(zhì)中地震波和電磁場的相互耦合以及震電波場在介質(zhì)中的傳播（具體公式詳見附錄A）.應(yīng)該指出，該理論假設(shè)雙相介質(zhì)中的孔隙是連通的并且被流體充滿，而實際孔隙介質(zhì)可能并非完全如此，所以Pride的震電控制方程組有其局限性.盡管如此，目前該理論依然是震電效應(yīng)研究中最受關(guān)注的一種理論.

基于Pride提出的震電波場宏觀控制方程組，Ren等［23］結(jié)合前人的研究成果，利用 Ge ＆ Chen［11］發(fā)展的直接計算廣義反透射系數(shù)的方法，得到了一套流體飽和層狀孔隙介質(zhì)中點源激發(fā)震電波場的數(shù)值模擬方法，并且基于點源疊加原理發(fā)展了有限斷層面源激發(fā)的震電波場的數(shù)值模擬方法.該方法首先是在柱坐標(biāo)下用一套基函數(shù)分解震電控制方程中所有的矢量和標(biāo)量，得到兩套方程組，分別對應(yīng)于SHTE模式和PSVTM模式（在SHTE模式中，SH波產(chǎn)生與它運動方向相同的電流，形成TE模式的電磁場；在PSVTM模式中，PSV波產(chǎn)生的電流位于PSV波運動面內(nèi)，形成TM模式的電磁場），具體公式詳見附錄B，然后根據(jù)位移－電磁場－應(yīng)力的連續(xù)性邊界條件得到各個界面上廣義反透射系數(shù)的具體表達(dá)式，進(jìn)而得到震電波場在頻率域的解.最后通過離散波數(shù)法［26］得到時間域的震電波場的解.

3 結(jié) 果

3.1 雙力偶點源激發(fā)的震電波場

本節(jié)利用上述基本理論和數(shù)值模擬方法，對均勻全空間中雙力偶點源激發(fā)的震電波場進(jìn)行數(shù)值模擬.如圖1所示x－o－y平面位于水平位置，z軸垂直向下.圖1（a、b）為雙力偶點源，其地震矩Mzx＝Mxz＝1.12×1014N·m.圖1（c、d）為有限斷層面源，大小為8km×5km，長度方向與x軸方向平行；斷層傾角為60°，滑移角為0°，地震矩為M0＝1.12×1017N·m.圖1所示模型中震源時間函數(shù)均采用峰值頻率為1Hz的Ricker子波，其時間偏移量為6s，震源中心位置均為（0，0，10km），接收點位于（3km，4km，0.05km）.圖1（a、c）為均勻半空間模型，圖1（b、d）為雙層介質(zhì)模型，不同模型下的介質(zhì)參數(shù)選取詳見表1.

圖2是均勻全空間模型（圖1a的半空間模型去掉自由表面）在雙力偶點源激發(fā)下產(chǎn)生的地震電磁信號，圖中自上而下分別表示位移場、電場、磁場的x、y和z分量.可清晰地分辨出由Ricker子波激發(fā)產(chǎn)生的P波和S波擾動，到達(dá)時間大約為9.035s和11.438s.伴隨地震波到達(dá)的電場信號在三個分量上都能清晰的看到.從圖2中可以看到S波激發(fā)同震電場信號的能力要明顯弱于P波.事實上，在我們的數(shù)值模擬過程中發(fā)現(xiàn)，對于一些介質(zhì)，S波激發(fā)的同震電場信號相較P波而言弱到幾乎可以忽略.總體來說，對于大多數(shù)的孔隙介質(zhì)，P波激發(fā)同震電場信號的能力要明顯強于 S波［8－9，13，20，23］，但是Gao ＆ Hu［20］指出對于某些實際的沉積巖，S波激發(fā)同震電場信號的振幅有可能超過P波激發(fā)的同震電場信號.由于P波不產(chǎn)生同震磁場信號，故僅在S波到達(dá)的時間處有磁場信號.

表1 均勻?qū)訝罱橘|(zhì)的參數(shù)Table 1 Properties of the layered porous media

為了進(jìn)一步考察均勻全空間中雙力偶點源激發(fā)的地震電磁信號的特征，我們計算了圖2對應(yīng)波場的頻率域響應(yīng)（圖3）.對于采用的峰值頻率為1Hz的Ricker子波，其頻譜成分主要集中在0～3Hz.從圖中可以看出產(chǎn)生的地震電磁信號頻譜成分與子波一致，其頻率集中在0～3Hz，并且震電波場九個分量都具有相近的優(yōu)勢頻率.結(jié)合圖2和圖3，位移場信號和電磁信號在到達(dá)時間和頻率上都有較好的一致性和相關(guān)性，其中到時一致性的特征與前人的研究結(jié)果［8－9，21－23］相吻合；湯吉等［4］開展了汶川地震強余震的同震電磁信號的實際觀測，也發(fā)現(xiàn)地震波和電磁場在到時方面具有較好的一致性.

3.2 不同模型的震電波場

本節(jié)將對雙力偶點源激發(fā)的全空間（圖1a的半空間模型去掉自由表面）、半空間（圖1a）和兩層介質(zhì)模型（圖1b），以及有限斷層面源作用下的半空間（圖1c）和兩層介質(zhì)模型（圖1d）進(jìn)行對比，分析模型結(jié)構(gòu)對地震電磁信號的影響.

圖4展示了全空間、半空間、雙層介質(zhì)模型中雙力偶點源和有限斷層面源激發(fā)產(chǎn)生的震電波場的Ex分量.全空間模型中雙力偶點源激發(fā)的Ex分量可清晰地分辨出在伴隨P波以及S波的電場信號先后到達(dá)接收點，并且伴隨P波的電場信號強于伴隨S波到達(dá)的電場.半空間模型中雙力偶點源激發(fā)的震電波場中，由于一般情況下P波產(chǎn)生同震電場信號的能力強于S波激發(fā)同震電場信號的能力，在自由界面處的Sv－P轉(zhuǎn)換波加強了在大約t＝11.438s時刻到達(dá)的電場信號（圖4b）.

與雙力偶點源相比，有限斷層面源產(chǎn)生的地震電磁信號在波形（圖4）和頻譜特征（圖5）方面均發(fā)生明顯變化，并且持續(xù)時間顯著增加（圖4）.此外，無論是雙力偶點源還是有限斷層面源，兩層介質(zhì)模型產(chǎn)生的地震電磁信號的波形都比均勻半空間模型的更為復(fù)雜，持續(xù)時間也更長（圖4）.

圖3 均勻全空間雙力偶點源激發(fā)的（a）地震波、（b）電場和（c）磁場信號頻譜示意圖Fig.3 Frequency spectrum of（a）seismic waves，（b）electric field，and （c）magnetic field generated by a double couple point source in a homogenous full space

3.3 不同類型震源時間函數(shù)的影響

圖6表示了Ricker子波、梯形函數(shù)、指數(shù)函數(shù)、Herrmann函數(shù)對應(yīng)的震源時間函數(shù).4種子波激發(fā)類型在時間域呈現(xiàn)不同的模式，為了便于比較，所設(shè)計的上述4種時間函數(shù)的時間窗都為2s，所選Ricker子波的峰值頻率為1Hz，其頻譜成分主要集中在0～3Hz，而梯形函數(shù)和指數(shù)函數(shù)采用漢寧窗濾波，濾波器的頻帶范圍為0～3Hz.

圖7表示了均勻半空間模型中雙力偶點源激發(fā)下，分別采取圖6所示的Ricker子波、梯形函數(shù)、指數(shù)函數(shù)、Herrmann函數(shù)［27］4種不同震源時間函數(shù)產(chǎn)生的地震電磁響應(yīng)，模型均采取相同的地震矩Mxz＝Mzx＝1.12×1014N·m，震源位置均為（0，0，10km），接收點均位于（3km，4km，0.05km）.圖7的結(jié)果顯示：圖6所示的不同震源時間函數(shù)激發(fā)的地震電磁信號的波形具有比較相似的整體特征，但具體的波形細(xì)節(jié)和振幅也存在一定的差異，其中，Ricker子波產(chǎn)生的信號幅值最大，其對應(yīng)的震電波場中可以清晰地分辨出伴隨P波以及S波的電場信號，而其他函數(shù)所產(chǎn)生的地震電磁輻射波形幅度相對較小，并且梯形函數(shù)和Herrmann函數(shù)對應(yīng)的震電波場相對較復(fù)雜.

圖4 雙力偶點源（a）全空間、（b）半空間、（c）兩層介質(zhì)以及有限斷層面源（d）半空間、（e）兩層介質(zhì)模型下的電場信號Ex示意圖Fig.4 Exof（a）full space，（b）half space，（c）two－layer model with a double couple source，and Exof（d）half space，（e）two－layer model with a finite fault source

圖5 雙力偶點源（a）全空間、（b）半空間、（c）兩層介質(zhì)以及有限斷層面源（d）半空間、（e）兩層介質(zhì)模型下的電場信號Ex頻譜示意圖Fig.5 Frequency spectrum of Ex（a）full space，（b）half space，（c）two－layer model with a double couple source，and Exof（d）half space，（e）two－layer model with a finite fault source

為了考察信號的頻譜特征，作為一個例子，圖8給出了上述Ricker子波、梯形函數(shù)、指數(shù)函數(shù)和Herrmann函數(shù)震源時間函數(shù)所產(chǎn)生的地震電磁信號x分量的頻譜圖.結(jié)果表明，不同震源時間函數(shù)激發(fā)的震電波場在頻率域的響應(yīng)呈現(xiàn)相似的整體特征，但在頻譜細(xì)節(jié)方面存在一定差異，并且在相同的地震矩情況下，震電波場頻譜的幅值變化明顯，其中Ricker子波幅值最大，產(chǎn)生的電磁場能量最強.

3.4 Ricker子波峰值頻率的影響

本小節(jié)重點考慮不同頻率激發(fā)的地震電磁信號的時域和頻域響應(yīng).上一小節(jié)的結(jié)果表明Ricker子波也能從整體上反映震電波場的主要特征，而且Ricker子波便于對頻率進(jìn)行定量分析，因此，本文以Ricker子波為例來考察震電波場對源的頻率響應(yīng)特征.圖9給出了均勻全空間峰值頻率分別為1，2，3Hz的Ricker子波激發(fā)的同震電磁信號.為清晰地顯示圖形中較弱的信號，對位移場P波所在的9～10.5s（圖9a方框部分）和S波產(chǎn)生的電場信號所屬的11.5～13.5s（圖9b方框部分）分別進(jìn)行了10倍的局部放大.可以看到：在某一頻率范圍內(nèi)，地震波隨著頻率的增大而顯著增強；伴隨P波和S波的電磁信號變化規(guī)律與地震波相似，峰值頻率越大，幅值越大，能量越集中（圖9（b，c））；不同頻率產(chǎn)生的同震電磁信號在相位和到達(dá)時間上保持一致.

圖8 震源時間函數(shù)為Ricker子波、梯形函數(shù)、指數(shù)函數(shù)和Herrmann函數(shù)對應(yīng)的（a）地震ux，（b）電場Ex，（c）磁場Hx的頻譜示意圖Fig.8 Frequency spectrum of（a）ux，（b）Ex，and（c）Hxinduced by source time function Ricker wavelet，trapezoid function，exponent function and Herrmann function

4 討論與結(jié)論

本研究基于天然地震頻段，重點探討了雙力偶點源以及有限斷層面源在不同的模型中所產(chǎn)生的震電波場的基本特征及可能存在的差異.數(shù)值模擬得到的地震波和電磁場在波形和頻譜方面具有較好的一致性和相關(guān)性.雖然前人研究很少涉及頻譜的討論，本文得到的頻譜方面的相關(guān)認(rèn)識尚有待更多理論模擬和實測資料的檢驗，但本文得到的波形的到時一致性這一特征與前人的研究結(jié)果相吻合，實際的野外同震電磁觀測結(jié)果也發(fā)現(xiàn)波場具有到時一致性的特征.全空間中伴隨P波以及S波都存在電場信號，但是S波激發(fā)同震電場信號的能力在一般情況下遠(yuǎn)弱于P波激發(fā)同震電場信號的能力.在半空間中由于自由界面的存在，Sv－P轉(zhuǎn)換波激發(fā)的同震電場信號與伴隨S波到達(dá)的同震電場信號疊加在一起，其幅值相比全空間模型顯著增強.

不同模型結(jié)構(gòu)的研究結(jié)果顯示，兩層介質(zhì)模型由于自由表面和介質(zhì)層間分界面的存在，產(chǎn)生了伴隨反射波和多次波傳播的電磁信號，故其波形較半空間模型更為復(fù)雜，持續(xù)時間更長，同時由于地震波在界面處發(fā)生反射，使得能量衰減，導(dǎo)致最終到達(dá)接收點的地震電磁信號有所減弱.總體而言，地震電磁信號的波形、振幅和持續(xù)時間等都不同程度受到模型結(jié)構(gòu)的影響，結(jié)構(gòu)越復(fù)雜，波形也越復(fù)雜.

對于有限斷層面源和雙力偶點源，前者可以視為斷層面內(nèi)的許多雙力偶點源的疊加，有限斷層面源的地震矩能量是隨著破裂傳播過程緩慢釋放的，而雙力偶點源的地震矩能量是同時釋放的，所以有限斷層面源產(chǎn)生的地震電磁信號的波型更復(fù)雜，持續(xù)時間也更長.

圖9 均勻全空間峰值頻率分別為1Hz、2Hz和3Hz的Ricker子波激發(fā)的（a）地震ux，（b）電場Ex，（c）磁場Hx示意圖（圖示中方框部分放大了10倍）Fig.9 （a）ux，（b）Ex，（c）Hxdue to the models with Ricker wavelet of 1Hz，2Hz，and 3Hz peak frequencies in a homogenous full space.The relative weak signals in the solid square have been amplified by a factor of 10for clear presentation

為了考察不同類型震源時間函數(shù)對地震和電磁波場的影響，我們分別選擇了Ricker子波、梯形函數(shù)、指數(shù)函數(shù)和Herrmann函數(shù)4種模型進(jìn)行了數(shù)值模擬，結(jié)果表明這4種不同模型激發(fā)的震電波場在波形和頻譜響應(yīng)方面整體上呈現(xiàn)相似的特征，但在波形和頻譜的細(xì)節(jié)方面存在一定差異，其中，Ricker子波產(chǎn)生的信號幅值最大，波形和頻譜結(jié)構(gòu)相對簡單、清晰.需要說明的是，雖然Ricker子波與地震的實際震源時間函數(shù)存在一定的差異，但其得到的結(jié)果既相對簡單、清晰，又能反映波場的整體特征而不失一般性，并且可以方便地考察結(jié)果的頻率響應(yīng)，所以，本文以Ricker子波為例對波場的頻率響應(yīng)進(jìn)行了探討，獲得了以下初步認(rèn)識，不同的Ricker子波峰值頻率會導(dǎo)致震電波場幅值的差異，但對信號的相位和到達(dá)時間并不產(chǎn)生顯著的影響.當(dāng)然，由于該認(rèn)識并無前人結(jié)果可以借鑒，其適用性尚有待更多數(shù)值模擬和實測資料的檢驗.

為了能更好地突出不同模型參數(shù)對震電波場的影響，本文選擇相對簡單的模型，分別對全空間、半空間和兩層介質(zhì)等模型結(jié)構(gòu)以及不同激發(fā)源（震源時間函數(shù)類型、Ricker子波峰值頻率）模型進(jìn)行了研究，獲得了震電波場的初步認(rèn)識.考慮到實際地球介質(zhì)結(jié)構(gòu)和震源模型可能比本研究采納的簡化模型更為復(fù)雜，今后有必要選擇更接近實際的模型進(jìn)一步進(jìn)行系統(tǒng)的數(shù)值模擬研究，為實際震電觀測資料的解釋提供理論依據(jù).

附錄A

對于均勻孔隙介質(zhì)，頻率域的震電控制方程組可表示如下（參見Pride［7］）：該方程組本質(zhì)上是Biot孔隙介質(zhì)理論和Maxwell方程組的結(jié)合，其中E和H 分別為電場和磁場，u和w分別為固體位移和流體滲濾位移，τ表示應(yīng)力張量，P表示流體壓力，F(xiàn)和f分別為介質(zhì)體力和作用在流體相上的壓力，ρ和ρf分別代表介質(zhì)平均密度和孔隙流體密度，L為隨頻率變化的動電耦合系數(shù)，J表示電流密度，C和M 分別表示外加的電流源和磁流源，KG為Gassmann體積模量，G為固體骨架剪切模量，KG與固體體積模量Ks、流體體積模量Kf和排除流體后骨架的體積模量Kfr具有如下的關(guān)系：

其中φ為介質(zhì)的孔隙度.

附錄B

對于底層為半無限空間的N＋1層層狀介質(zhì)模型，第j（1≤j≤N）層的邊界記為z＝z（j－1）和z＝z（j），頂層可以是半空間，也可以是一厚度有限的層（此時將其自由表面記為z＝z（0）），假設(shè)源位于第s層并且2≤s≤N，即源不在最頂層也不在最底層.

在z軸垂直向下的柱坐標(biāo)下，用一套基函數(shù)（見Ren等［23］的附錄A）分解震電控制方程中所有的矢量和標(biāo)量，得到兩套方程組，分別對應(yīng)于SHTE模式和PSVTM模式，他們都可以寫成如下的形式：

其中，A（j）代表第j層介質(zhì)物性參數(shù)所控制的系數(shù)矩陣，F(xiàn)為外加的力源項（其具體表達(dá)式可參考Ren等［23］的附錄B），對于SHTE模式和PSVTM模式，分別有

同時對于電磁場垂直分量，有

（References）

［1］Matsushima M， Honkura Y， Oshiman N， et al.Seismoelectromagnetic effect associated with the I·zmit earthquake and its aftershocks.Bull.Seism.Soc.Am.，2002，92：350－360.

［2］Karakelian D，Beroza G C，Klemperer S L，et al.Analysis of ultralow－frequency electromagnetic field measurements associated with the 1999 M7.1Hector Mine，California，earthquake sequence.Bull.Seism.Soc.Am.，2002，92：1513－1524.

［3］湯吉，詹艷，王立鳳等.5月12日汶川8.0級地震強余震觀測的電磁同震效應(yīng).地震地質(zhì)，2008，30（3）：739－745.Tang J，Zhan Y，Wang L F，et al.Coseismic signal associated with aftershock of the Ms8.0Wenchuan earthquake.Seismology and Geology （in Chinese），2008，30（3）：739－745.

［4］湯吉，詹艷，王立鳳等.汶川地震強余震的電磁同震效應(yīng).地球物理學(xué)報，2010，53（3）：526－534.Tang J，Zhan Y，Wang L F，et al.Electromagnetic coseismic effect associated with aftershock of Wenchuan Ms8.0 earthquake.Chinese J.Geophys.（in Chinese），2010，53（3）：526－534.

［5］Ivanov A G.Effect of electrization of earth layers by elastic waves passing through them.Dokl.Akad.Nauk.SSSR，1939，24：42－45.

［6］Frenkel J.On the theory of seismic and seismoelectric phenomena in a moist soil.J.Eng.Mech.，1944，8（4）：230－241.

［7］Pride S R.Governing equations for the coupled electromagnetics and acoustics of porous media.Phys.Rev.B，1994，50：15678－15696.

［8］Haartsen M W，Pride S R.Electroseismic waves from point sources in layered media.J.Geophys.Res.，1997，102：24745－24769.

［9］Garambois S， Dietrich M. Full waveform numerical simulations of seismoelectromagnetic wave conversions in fluid－saturated stratified porous media.J.Geophys.Res.，2002，107：40－58.

［10］Chen X F.Seismogram synthesis in multi－layered half－space part I.Theoretical formulations.Earthquake Research in China，1999，13：149－174.

［11］Ge Z X，Chen X F.Wave propagation in irregularly layered elastic models：A boundary element approach with a global reflection／transmission matrix propagator.Bull.Seism.Soc.Am.，2007，97（3）：1025－1031.

［12］Ren H X，Huang Q H，Chen X F.A new numerical technique for simulating the coupled seismic and electromagnetic waves in layered porous media.Earthq.Sci.，2010，23：167－176.

［13］Ren H，Huang Q H，Chen X.Analytical regularization of the high－frequency instability problem in numerical simulation of seismoelectric wave－fields in multi－layered porous media.Chinese J.Geophys.，2010（b），53：506－511.

［14］Martner S T，Sparks N R.The electroseismic effect.Geophys.，1959，24：297－308.

［15］Thompson A H，Gist G A.Geophysical applications of electrokinetic conversion.Leading Edge，1993，12：1169－1173.

［16］Thompson A H，Hornbostel S，Burns J，et al.Field tests of electroseismic hydrocarbon detection.Geophys.，2007，72：N1－N9.

［17］Zhu Z，Cheng C H，Toks?z M N.Electro－kinetic conversion in a fluid－saturated porous rock sample.SEG Technical Program Expanded Abstracts，1994，1057－1060.

［18］Zhu Z，Haartsen M W，Toks?z M N.Experimental studies of electrokinetic conversions in fluid－saturated borehole models.Geophys.，1999，64：1349－1356.

［19］Zhu Z，Haartsen M W，Toks?z M N.Experimental studies of seismoelectric conversions in fluid－saturated porous media.J.Geophys.Res.，2000，105：28055－28064.

［20］Huang Q H. Rethinking earthquake－related DC－ULF electromagnetic phenomena：towards a physics－based approach.Nat.Haiards Earth Syse.Sci.，2011，11（11）：2941－2949.

［21］Gao Y H，Hu H S.Seismoelectromagnetic waves radiated by a double couple source in a saturated medium.Geophys.J.Int.，2010，181：873－896.

［22］Hu H S，Gao Y X.Electromagnetic field generated by a finite fault due to electrokinetic effect.J.Geophys.Res.，2011，116.

［23］Ren H X，Chen X F，Huang Q H.Numerical simulation of coseismic electromagnetic fields associated with seismic waves due to finite faulting in porous media.Geophys.J.Int.，2012，188：925－944.

［24］Gao Y X，Chen X F，Hu H S，et al.Early electromagnetic waves from earthquake rupturing： I. theoretical formulations.Geophys.J.Int.，2013，192：1288－1307.

［25］Gao Y X，Chen X F，Hu H S，et al.Early electromagnetic waves from earthquake rupturing：II.validation and numerical experiments.Geophys.J.Int.，2013，192：1308－1323.

［26］Bouchon M.A simple method to calculate Green′s functions for elastic layered media.Bull.Seism.Soc.Am.，1981，71：959－971.

［27］Herrmann R B.SH－wave generation by dislocation source—a numerical study.Bull.Seism.Soc.Am.，1979，69：1－15.