三維正交各向異性孔隙介質震電波場數值模擬研究

穆納爾丁·托合提， 王一博， 肖文交,， 周可法

1 中國科學院新疆生態與地理研究所， 荒漠與綠洲生態國家重點實驗室， 烏魯木齊 830011 2 新疆礦產資源與數字地質重點實驗室， 烏魯木齊 830011 3 中國科學院新疆礦產資源研究中心， 烏魯木齊 830011 4 中國科學院大學地球與行星科學學院， 北京 100049 5 中國科學院地質與地球物理研究所, 北京 100029

0 引言

在含流體的孔隙介質中地震波和電磁波可以互相轉換，即地震波在含流體的孔隙介質傳播過程中可以產生電磁波，同樣電磁波在含流體的孔隙介質傳播過程中可以激發地震波，這種耦合現象叫震電效應.震電信號的產生與介質物性參數密切相關，研究震電信號有助于了解震電效應的產生機理以及傳播過程，具有重要的理論意義和實際價值.震電效應的研究最早可以追溯到1939年，蘇聯的Ivanov發現，在多孔介質中地震波可以引起電勢變化(Ivanov，1939).隨后Frenkel(1944)用動電理論解釋了這種現象.后來Biot(1955，1956a，b)在Frenkel理論基礎之上發展了一套孔隙介質中描述地震波傳播的波動方程理論，為后來Pride震電耦合方程組的建立奠定了基礎(Pride，1994).

目前有三種耦合的震電信號已經成為研究的熱點(Wang et al.，2020)，分別是伴隨震電信號(Haartsen and Pride，1997； Gao and Hu，2010)、界面震電信號(Haartsen and Pride，1997； Garambois and Dietrich，2002；Haines and Pride，2006；Thompson et al.，2005；Grobbe and Slob，2016)和隱式震電信號(Ren et al.，2016a).作為局部響應，伴隨震電信號是由相對的固液運動引起的，并以與地震波相同的速度傳播.當地震波遇到兩種不同的多孔介質之間的界面或者多孔介質與另一種介質(例如流體)的界面時，會產生界面震電信號(Mikhailov et al.，1997),并以高于地震波的速度傳播.Ren等(2016a)發現，地震波在含流體的層狀孔隙介質中，除了產生伴隨震電信號和界面震電信號之外，還會產生一種新的電磁波叫隱失電磁波(evanescent electromagnetic waves),當地震波傳到界面且入射角大于臨界角時，在界面處會產生隱式電磁波.隱式震電信號其實是當地震波以超過臨界角入射產生的非均勻電磁波，本質上仍是界面震電信號，只是該非均勻電磁波不能向全空間任意方向傳播了，只能沿著界面以相應地震波的速度傳播，表現類似“同震信號”.Martner和Sparks(1959)通過野外實驗研究了放在不同深度的爆炸源產生的震電信號，結果表明在地表可以觀測到有一定埋深的分化層底部的界面產生的震電信號.嚴洪瑞等(1999)在我國的大慶油田做了野外試驗，進行了震電勘探，成功地觀測到了產生的震電信號，但是發現產生的震電信號信噪比較低，信號也比較弱.考慮到震電轉換信號本身較弱且隨著距離的增加呈指數衰減，有學者提出將接收點置于測量目標區域附近，即井中震電勘探的想法(陳本池等，2003；張元中等，2005).還有很多學者也進行了一系列野外和室內實驗，主要驗證了界面電磁波和伴隨震電信號的存在(胡恒山等，2001；Garambois and Dietrich，2002；陳本池，2007；Bordes et al.，2006； Zhu et al.，2008，2016；Zhu and Toks?z, 2013; Schakel et al.，2011a，b，2012；Schoemaker et al.，2012； Bordes et al.，2006，2015；Roubinet et al.，2016；Butler et al.，2018；Peng et al.，2017，2019).

為了定量描述不同類型的震電信號的傳播特征, Pride(1994)利用動電理論推導出了一組震電控制方程.控制方程組由兩部分組成，分別是描述電磁波傳播和擴散的麥克斯韋方程組和描述孔隙介質中彈性波傳播的Biot方程組(Biot，1955，1956a，b).基于Pride模型，很多學者進行了一系列數值模擬和理論研究，探討了不同介質中產生的震電波場特征，包括均勻介質(Pride and Haartsen，1996；Gao and Hu，2010；Ren et al.，2012；Slob and Mulder， 2016)、水平層狀介質(Haartsen and Pride，1997；Garambois and Dietrich，2002；Jardani et al.，2010；Gao and Hu，2010；Gao et al., 2013a,b，2017a,b，2019; Jardani and Revil，2015；Ren et al.，2010，2016a，b；Zyserman et al.，2015；Grobbe and Slob，2016；Jouniaux and Zyserman，2016；Monachesi et al.，2018a，b)和井孔模型(胡恒山和王克協，1999；Guan and Hu，2008；Guan et al.，2013).Haartsen和Pride(1997)用數值方法確定了層狀孔隙介質中點源的電震響應.Garambois和Dietrich(2002)基于Pride (1994)耦合方程給出了在層狀飽和孔隙介質中震電耦合波場的模擬技術.Hu和Liu(2002)通過準靜態近似方法，模擬了聲電效應測井響應.Gao和Hu(2010)研究了剪切源在層狀飽和孔隙介質中的震電波場響應.Haines和Pride(2006)基于有限差分算法，計算了準靜態近似下層狀孔隙介質的震電波場特征.考慮到震電波場中準靜態近似會帶來誤差，Gao等(2017b)的研究工作對誤差分析做了相應的探討.Guan等(2018)基于準靜態近似的數值模擬研究，分析了孔彈性波傳播特征對震電波場的影響.Revil和Mahardika(2013)提出了另一震電耦合方程組，描述了在不飽和多孔介質中震電波場的傳播，其與Pride的區別在于動電耦合系數L的表達式不同，耦合系數L依賴于凈剩離子的含量.基于Revil和Mahardika(2013)的震電耦合模型，Save等人進行了數值模擬，探討了震電耦合效應(Sava et al.，2014).

通過回顧國內外相關的研究進展，可以看到目前關于孔隙介質中的震電波場數值模擬研究主要是基于各向同性介質展開的，Tohti等(2020)之前做過有關各向異性孔隙介質中的震電波場數值模擬研究，但是只考慮了二維垂直對稱軸的橫向各向同性類型.為了更為全面地認識實際孔隙介質中產生的震電效應特征，需要開展三維正交各向異性孔隙模型的數值模擬研究.

本文是對我們之前工作的擴展(Tohti et al.， 2020)，研究了三維正交各向異性介質的震電耦合方程組，通過設計幾種不同的正交各向異性模型，模擬震電信號在正交各向異性介質中的傳播特征，分析了介質參數對震電波場的影響.

1 方法原理

1.1 各向異性介質震電耦合方程

基于Pride的控制方程(Pride, 1994)，Slob和Mulder在頻率-波數域給出了三組基本方程(Slob and Mulder, 2016).它們分別是修正的麥克斯韋方程、運動方程和本構方程.在準靜態近似下，我們可以將運動方程寫成如下形式：

(1)

(2)

本構方程在時間域可以寫成如下形式：

(3)

(4)

Slob和Mulder(2016)給出了如方程(5)所示的麥克斯韋方程組：

(5)

(6)

根據前人的研究結果 (Biot, 1956a,b; Carcione, 1996), 方程(1)—(6)在正交各向異性孔隙介質中可以寫成如下形式：

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

1.2 震電耦合方程組有限差分離散格式

本文用有限差分算法對震電耦合方程組進行離散(Virieux, 1984; Tohti et al.，2020).用海綿吸收邊界條件吸收邊界上的反射波場(Petropoulos et al., 1998).根據以下原理，對方程(7)—(22)中的偏導數進行近似，

(23)

根據方程 (23), 可得到耦合方程組的有限差分離散格式如下:

(25)

+C1(x1,x2,x3)[w2(g3)-w2(g4)]+C1(x1,x2,x3)[w3(g5)-w3(g6)]]

(30)

+C2(x1,x2,x3)[w2(g3)-w2(g4)]+C2(x1,x2,x3)[w3(g5)-w3(g6)]],

(31)

+C3(x1,x2,x3)[w2(g3)-w2(g4)]+C3(x1,x2,x3)[w3(g5)-w3(g6)]],

(32)

×[v3(g3)-v3(g4)]],

(33)

×[v3(g1)-v3(g2)]],

(34)

×[v1(g3)-v1(g4)]],

(35)

×[v2(g3)-v2(g4)]+C3(x1,x2,x3)[v3(g5)-v3(g6)]+M(x1,x2,x3)[w1(g1)-w1(g2)]

+M(x1,x2,x3)[w2(g3)-w2(g4)]+M(x1,x2,x3)[w3(g5)-w3(g6)]],

(36)

(37)

(38)

(39)

2 數值模擬

2.1 方法驗證

為了驗證我們數值模擬算法的正確性，第一步，將各向異性參數設置為ε1=ε2=0,δ1=δ2=0,γ1=γ2=0,δ3=0來模擬各向同性介質.

本文使用Slob和Mulder(2016)的3D解析解來驗證了模擬代碼，使用1 kHz Ricker子波作為震源時間函數，模型大小為624×624×624網格數，所有方向的網格間距均為0.1 m，時間步長為0.02 ms.我們將爆炸源置于模型的正中心，并在接收器R處接收信號.源和接收器的坐標如圖1所示.

在所有模擬例子中，我們都對均勻介質進行了模擬.使用表1中的模型A的參數測試了代碼.數值解與解析解的對比結果如圖2所示.從對比結果可以看出，FDTD數值解與解析解基本吻合，這說明本文的模擬算法是正確的.

圖1 接收器和源的坐標示意圖Fig.1 Configuration of receiver and source location

圖2 震電解析解和數值解之間對比 (a), (b), (c) 分別代表x1, x2和x3方向上的固體速度分量(V)； (d), (e), (f) 跟 (a), (b), (c)一樣，但表示相對速度分量(W)； (g), (h), (i) 跟 (a), (b), (c)一樣，但表示電場分量(E).Fig.2 Comparison between analytical solutions and numerical solutions of seismoelectric waves (a), (b), (c) represent the solid velocity (V) components on x1, x1, and x1 directions, respectively; (d), (e), (f) are the same as (a), (b), (c), respectively, but for filtration velocity (W) components; (g), (h), (i) are the same as (a), (b), (c), respectively, but for electric field (E) components.

2.2 數值例子

在本節中，我們分別模擬爆炸源在模型A、模型B、模型C、模型D、模型E和模型F中產生的震電響應.其中模型B、模型C和模型D只考慮了骨架的各向異性來分析骨架各向異性對震電波場的影響.模型F考慮彎曲度和滲透率的各向異性.我們將主頻為1 kHz 的Ricker子波作為震源時間函數.模型尺寸為624×624×624網格數，空間步長為Δx=0.1 m，時間步長為0.02 ms.我們在源上方x2方向上布置了100個接收器，接收器與源位置之間的最小距離為15 m.

首先用各向同性模型A來模擬震電波場，模型A的參數如表1所示.圖3中顯示了在10 ms時刻的波場快照.我們觀察到爆炸源在均勻各向同性多孔介質中激發兩種類型的地震波.它們分別是快縱波(Pf)和慢縱波(PS).這兩種波都產生伴隨震電場.快速縱波類似于在彈性介質中產生的縱波.然而，慢縱波的傳播特征取決于流體性質.慢縱波在流體黏度的影響下傳播的慢，以靜態模式出現在震源位置，如圖3所示，其中EPf和EPs分別是Pf和Ps的伴隨震電場.

為了研究震電信號在正交各向異性多孔介質中的傳播特性，并分析流體黏度、固體骨架各向異性和彎曲度對震電波傳播的影響，我們另外設置了具有不同參數的5種模型(標記為B到F)，如表1所示.

表1 在數值模擬中使用的模型參數Table 1 Model parameters used in numerical modeling

模型B與模型A的區別在于，模型A是各向同性的，而模型B是正交各向異性介質.我們在圖4中顯示了在10 ms時刻的模擬結果.可以觀察到，除了快速波和慢速縱波外，爆炸源在各向異性孔隙介質中還會產生橫波(S)(Ben-Menahem and Sena, 1990；Ben-Menahem et al., 1991).我們發現，波的各向異性在不同的平面上是不同的.存在[x1,x2]平面中的快S波和慢S波.ESf和ESS分別代表快和慢S波的伴隨震電場.

為了分析固體骨架各向異性對震電信號傳播的影響，我們模擬了模型C的震電響應.與模型B相比，模型C中各向異性參數ε2的值增加了，可以分析參數ε2對震電波場傳播的影響.根據正交各向異性介質的定義，我們知道ε2控制[x1,x3]平面中波場的各向異性，而ε1控制[x2,x3]平面中波場的各向異性.因此我們只需要關注[x1,x3]平面中波場的變化就能分析出ε2對波場傳播的影響.在10 ms時刻的模擬結果如圖5所示.從圖5中可以看到，在模型C產生的 [x1,x3]平面中的Pf波比模型B在 [x1,x3]平面中的Pf波顯示出更大的各向異性效應.這說明各向異性參數ε(ε1,ε2)的值主要控制Pf波的各向異性(Tsvankin,1997; Wang et al., 2018).由于在模型C中ε2的值大于ε1的值，因此在模型C產生的[x1,x3]平面中的Pf波比模型C在[x2,x3]平面里的Pf波表現出更大的各向異性.

為了分析δ對震電信號傳播的影響，我們模擬了模型D的震電響應.模型D與模型C的區別在于，模型D中的δ2值發生了變化.在正交各向異性介質中，我們知道δ2和δ1分別控制[x1,x3]和[x2,x3]平面里波場的各向異性，而δ3控制[x1,x2]平面里波場的各向異性.因此我們只需要關注[x1,x3]平面中波場的變化就能分析出δ2對波場傳播的影響.圖6顯示了10 ms時刻的波場快照.我們觀察到兩種橫波(S)在[x1,x3]平面中彼此分離了，換句話說，在[x1,x3]平面中的S波是在δ2的作用下分裂了，而在[x2,x3]和[x1,x2]平面中橫波的變化不大，因為δ1和δ3的值沒有發生變化.由此可以確定δ(δ1,δ2,δ3)對S波的各向異性的影響要大于對Pf波的影響(Wang et al., 2018).由于δ1,δ2和δ3的值不同，在模型B、模型C和模型D中產生的S波在不同平面上表現出不同的各向異性.

為了分析流體黏度對震電信號傳播的影響，我們模擬了模型E的震電響應.模型E與模型D的區別在于，模型E中的流體黏度值減小了.在10 ms時刻的模擬結果如圖7所示.我們觀察到慢縱波(PS)及其伴隨震電波(EPs)開始在介質中傳播，換句話說，當流體黏度值減少時，慢縱波開始以傳播模式出現了.由此證實了流體黏度對慢縱波及其伴隨震電場的傳播有較大的影響，當流體黏度降低時，慢縱波的傳播速度增加.

最后，我們模擬了模型F的震電響應，以分析彎曲度對震電波傳播的影響.模型F和模型E的區別在于，在模型F中彎曲度是各向異性的.由于耦合系數與彎曲度有關，耦合系數也是各向異性的.從Guo(2012)中可以看出，彎曲度和滲透率之間存在相關性，因此模型F中的滲透率也是各向異性的.圖8顯示了10 ms時刻的模擬結果.我們觀察到，慢縱波及其伴隨震電波(EPs)的波場快照形狀顯示為橢圓形，尤其是在[x1,x3]平面中，而模型E的模擬結果中慢縱波的波場快照顯示為圓形.由此可以確定，彎曲度對慢縱波及其伴隨震電場的傳播有影響.慢縱波在[x1,x3]平面中具有更大的各向異性，這是因為x1和x3方向之間的彎曲度差異大于其他方向之間的差異.

2.3 運算效率

本文用Open Mp將FDTD代碼并行化，在15 ms的模擬時間下用不同網格數和線程數測試了并行計算的性能.表2中給出了模型大小和線程數.基于表2，可以分析加速比和計算時間隨著線程數的變化.

圖9顯示了不同模型大小和線程數下的加速比和運行時間.我們觀察到，隨著線程數量的增加，運行時間減少，加速比增加.

表2 不同模型大小和線程數下的運行時間Table 2 Run times for different model sizes and different number of threads

圖3 模型A的震電波場快照Fig.3 Seismoelectric snapshot of Model A

圖4 模型B的震電波場快照Fig.4 Seismoelectric snapshot of Model B

圖5 模型C的震電波場快照Fig.5 Seismoelectric snapshot of Model C

圖6 模型D的震電波場快照Fig.6 Seismoelectric snapshot of Model D

圖7 模型E的震電波場快照Fig.7 Seismoelectric snapshot of Model E

圖8 模型F的震電波場快照Fig.8 Seismoelectric snapshot of Model F

圖9 在不同模型大小和不同線程數下的運行時間(a)和加速比(b)Fig.9 Variation of run times(a) and speedup (b) with the number of threads for different model sizes

3 結論和認識

本文研究了各向異性孔隙介質震電耦合方程組，推導了三維正交各向異性孔隙介質情況下的時間域交錯網格有限差分離散格式，模擬了爆炸源的震電響應，分析了介質參數對震電信號的影響，得出了如下結論.

(1) 爆炸源在正交各向異性的多孔介質中產生四種波，分別是快縱波(Pf)，慢縱波(PS)和兩個可分離的橫波.這四種地震波在均勻的多孔介質中都能產生伴隨震電場.流體黏度和彎曲度對慢縱波的傳播有影響.在流體黏度的影響下，慢縱波緩慢傳播.彎曲度各向異性時，慢縱波的波場快照形狀呈橢圓形，而彎曲度各向同性時，慢縱波的波場快照形狀呈圓形.Pf波和S波的各向異性分別對ε(ε1,ε2)和δ(δ1,δ2,δ3)的值更敏感.

(2) 本文使用海綿吸收邊界條件來壓制了人工邊界上的反射波，邊界條件有效地吸收了邊界處的反射信號.

(3) 基于Open Mp將模擬代碼并行化，發現并行效果明顯使運算速度提高了.

致謝感謝編輯和兩位審稿人的寶貴意見.