三維正交各向異性孔隙介質(zhì)震電波場數(shù)值模擬研究

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

1 中國科學(xué)院新疆生態(tài)與地理研究所， 荒漠與綠洲生態(tài)國家重點實驗室， 烏魯木齊 830011 2 新疆礦產(chǎn)資源與數(shù)字地質(zhì)重點實驗室， 烏魯木齊 830011 3 中國科學(xué)院新疆礦產(chǎn)資源研究中心， 烏魯木齊 830011 4 中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院， 北京 100049 5 中國科學(xué)院地質(zhì)與地球物理研究所, 北京 100029

0 引言

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

目前有三種耦合的震電信號已經(jīng)成為研究的熱點(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).作為局部響應(yīng)，伴隨震電信號是由相對的固液運動引起的，并以與地震波相同的速度傳播.當(dāng)?shù)卣鸩ㄓ龅絻煞N不同的多孔介質(zhì)之間的界面或者多孔介質(zhì)與另一種介質(zhì)(例如流體)的界面時，會產(chǎn)生界面震電信號(Mikhailov et al.，1997),并以高于地震波的速度傳播.Ren等(2016a)發(fā)現(xiàn)，地震波在含流體的層狀孔隙介質(zhì)中，除了產(chǎn)生伴隨震電信號和界面震電信號之外，還會產(chǎn)生一種新的電磁波叫隱失電磁波(evanescent electromagnetic waves),當(dāng)?shù)卣鸩▊鞯浇缑媲胰肷浣谴笥谂R界角時，在界面處會產(chǎn)生隱式電磁波.隱式震電信號其實是當(dāng)?shù)卣鸩ㄒ猿^臨界角入射產(chǎn)生的非均勻電磁波，本質(zhì)上仍是界面震電信號，只是該非均勻電磁波不能向全空間任意方向傳播了，只能沿著界面以相應(yīng)地震波的速度傳播，表現(xiàn)類似“同震信號”.Martner和Sparks(1959)通過野外實驗研究了放在不同深度的爆炸源產(chǎn)生的震電信號，結(jié)果表明在地表可以觀測到有一定埋深的分化層底部的界面產(chǎn)生的震電信號.嚴(yán)洪瑞等(1999)在我國的大慶油田做了野外試驗，進(jìn)行了震電勘探，成功地觀測到了產(chǎn)生的震電信號，但是發(fā)現(xiàn)產(chǎn)生的震電信號信噪比較低，信號也比較弱.考慮到震電轉(zhuǎn)換信號本身較弱且隨著距離的增加呈指數(shù)衰減，有學(xué)者提出將接收點置于測量目標(biāo)區(qū)域附近，即井中震電勘探的想法(陳本池等，2003；張元中等，2005).還有很多學(xué)者也進(jìn)行了一系列野外和室內(nèi)實驗，主要驗證了界面電磁波和伴隨震電信號的存在(胡恒山等，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)利用動電理論推導(dǎo)出了一組震電控制方程.控制方程組由兩部分組成，分別是描述電磁波傳播和擴散的麥克斯韋方程組和描述孔隙介質(zhì)中彈性波傳播的Biot方程組(Biot，1955，1956a，b).基于Pride模型，很多學(xué)者進(jìn)行了一系列數(shù)值模擬和理論研究，探討了不同介質(zhì)中產(chǎn)生的震電波場特征，包括均勻介質(zhì)(Pride and Haartsen，1996；Gao and Hu，2010；Ren et al.，2012；Slob and Mulder， 2016)、水平層狀介質(zhì)(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)和井孔模型(胡恒山和王克協(xié)，1999；Guan and Hu，2008；Guan et al.，2013).Haartsen和Pride(1997)用數(shù)值方法確定了層狀孔隙介質(zhì)中點源的電震響應(yīng).Garambois和Dietrich(2002)基于Pride (1994)耦合方程給出了在層狀飽和孔隙介質(zhì)中震電耦合波場的模擬技術(shù).Hu和Liu(2002)通過準(zhǔn)靜態(tài)近似方法，模擬了聲電效應(yīng)測井響應(yīng).Gao和Hu(2010)研究了剪切源在層狀飽和孔隙介質(zhì)中的震電波場響應(yīng).Haines和Pride(2006)基于有限差分算法，計算了準(zhǔn)靜態(tài)近似下層狀孔隙介質(zhì)的震電波場特征.考慮到震電波場中準(zhǔn)靜態(tài)近似會帶來誤差，Gao等(2017b)的研究工作對誤差分析做了相應(yīng)的探討.Guan等(2018)基于準(zhǔn)靜態(tài)近似的數(shù)值模擬研究，分析了孔彈性波傳播特征對震電波場的影響.Revil和Mahardika(2013)提出了另一震電耦合方程組，描述了在不飽和多孔介質(zhì)中震電波場的傳播，其與Pride的區(qū)別在于動電耦合系數(shù)L的表達(dá)式不同，耦合系數(shù)L依賴于凈剩離子的含量.基于Revil和Mahardika(2013)的震電耦合模型，Save等人進(jìn)行了數(shù)值模擬，探討了震電耦合效應(yīng)(Sava et al.，2014).

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

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

1 方法原理

1.1 各向異性介質(zhì)震電耦合方程

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

(1)

(2)

本構(gòu)方程在時間域可以寫成如下形式：

(3)

(4)

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

(5)

(6)

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

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

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

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

(23)

根據(jù)方程 (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 數(shù)值模擬

2.1 方法驗證

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

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

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

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

圖2 震電解析解和數(shù)值解之間對比 (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 數(shù)值例子

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

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

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

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

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

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

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

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

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

2.3 運算效率

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

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

表2 不同模型大小和線程數(shù)下的運行時間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 在不同模型大小和不同線程數(shù)下的運行時間(a)和加速比(b)Fig.9 Variation of run times(a) and speedup (b) with the number of threads for different model sizes

3 結(jié)論和認(rèn)識

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

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

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

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

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