基于Low rank近似的黏彈性介質(zhì)衰減補(bǔ)償微地震逆時(shí)定位與裂縫成像

唐杰, 劉英昌， 李聰， 孫成禹

中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 山東青島 266580

0 引言

微地震監(jiān)測(cè)中震源定位是其核心問(wèn)題(Maxwell and Urbancic，2001)，在彈性介質(zhì)中通過(guò)波場(chǎng)逆時(shí)反傳微地震信號(hào)能夠有效確定震源位置(Steiner et al.，2008；Larmat et al.，2009)，而對(duì)于黏彈性介質(zhì)如果采用常規(guī)逆時(shí)反傳，由于地下介質(zhì)具有黏彈性會(huì)發(fā)生吸收衰減作用，引起地震波振幅減弱、相位失真，無(wú)法確定震源真實(shí)位置，因此需要對(duì)波場(chǎng)進(jìn)行衰減補(bǔ)償，采用逆時(shí)定位方法獲取真實(shí)的震源位置.

黏彈性介質(zhì)中微地震震源逆時(shí)定位包括兩個(gè)步驟：在時(shí)間上反傳檢波器接收到的地震數(shù)據(jù)與補(bǔ)償黏彈性介質(zhì)對(duì)地震波的吸收衰減作用.黏彈性介質(zhì)的速度頻散和能量耗散可以采用品質(zhì)因子Q描述，Kjartansson(1979)提出了常Q模型，以參考速度、衰減因子和參考頻率為參數(shù)，衰減因子是頻率的線性函數(shù).進(jìn)行波場(chǎng)逆時(shí)反傳的理論基礎(chǔ)是地震波正演模擬，近年來(lái)眾多學(xué)者對(duì)黏聲/黏彈性介質(zhì)正演模擬的算法進(jìn)行了研究.頻率域數(shù)值模擬基于單頻率對(duì)空間網(wǎng)格點(diǎn)求解，各頻率獨(dú)立計(jì)算不存在誤差累積，可以方便的進(jìn)行并行計(jì)算，但是在大型3D模型中該算法無(wú)法適應(yīng)消息傳遞接口(Message Passing Interface，MPI)并行算法及圖形處理器(Graphics Processing Unit，GPU)加速導(dǎo)致效率降低；時(shí)間域數(shù)值模擬耗費(fèi)內(nèi)存小，操作靈活，但是存在著時(shí)間循環(huán)累積的誤差.Carcione等(Carcione et al.，2002；Carcione，2009)通過(guò)GL(Grünwald-Letnikov)近似求解時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)，實(shí)現(xiàn)了黏聲及黏彈性介質(zhì)波動(dòng)方程數(shù)值模擬，但該方法在求解時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)時(shí)，需要保存每一時(shí)刻的應(yīng)力、應(yīng)變值，計(jì)算效率低，且內(nèi)存需求較高.楊仁虎等(2009)提出采用拉梅差異矩陣簡(jiǎn)化方程，通過(guò)該方法模擬黏彈性介質(zhì)波場(chǎng)具有比GL近似更高的效率.Carcione(2010)提出基于時(shí)間為二階導(dǎo)數(shù)、空間為分?jǐn)?shù)階的黏聲波動(dòng)方程，該方程將分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)轉(zhuǎn)化到分?jǐn)?shù)階空間導(dǎo)數(shù)上，有效減少了內(nèi)存占用，解決GL近似的效率問(wèn)題.黏聲方程中耗散項(xiàng)與頻散項(xiàng)的解耦促進(jìn)了衰減補(bǔ)償逆時(shí)延拓的發(fā)展，Treeby和Cox(2010)基于分?jǐn)?shù)階黏聲波動(dòng)方程，將耗散項(xiàng)與頻散項(xiàng)解耦，提高了計(jì)算效率，解耦的耗散項(xiàng)與頻散項(xiàng)也可以用于波場(chǎng)反傳中地震波的衰減補(bǔ)償計(jì)算，有利于震源定位和偏移成像.Zhu和Carcione(2014)利用空間分?jǐn)?shù)階導(dǎo)數(shù)將黏彈性介質(zhì)波動(dòng)方程耗散與頻散解耦.吳玉等(2017)將解耦的分?jǐn)?shù)階拉普拉斯算子用于黏聲介質(zhì)逆時(shí)偏移，改善了深層構(gòu)造的成像精度.對(duì)于波數(shù)-空間域的混合域算子求解時(shí)，計(jì)算量非常大，F(xiàn)omel等(2013)提出low rank 近似混合域算子，減少傅里葉變換次數(shù)提高計(jì)算效率.Song等(2013)結(jié)合low rank頻散低和有限差分法計(jì)算效率高的優(yōu)點(diǎn)，提出了low rank有限差分法解決聲波正演問(wèn)題.Sun等(2015)利用low rank近似模擬黏聲介質(zhì)波場(chǎng)，實(shí)現(xiàn)了黏聲介質(zhì)衰減補(bǔ)償逆時(shí)偏移成像；袁雨欣等(2018)采用low rank算法求解解耦的彈性波方程，有效地減少了波場(chǎng)模擬數(shù)值頻散，提高模擬精度.黃金強(qiáng)和李振春(2019)在各向異性介質(zhì)中應(yīng)用low rank近似進(jìn)行正演及逆時(shí)偏移，提高計(jì)算效率.近年來(lái)，為解決分?jǐn)?shù)階拉普拉斯算子計(jì)算復(fù)雜介質(zhì)時(shí)成本高的問(wèn)題，眾多學(xué)者不斷改進(jìn)算法進(jìn)行黏彈介質(zhì)的數(shù)值模擬(Chen et al.，2016；Yao et al.，2017；Wang et al.，2018；Guo et al.，2019).

在微地震震源成像條件的研究中，McMechan(1982)提出了最大振幅成像條件定位震源，但是由于需要不斷檢索波場(chǎng)最大值，該算法計(jì)算效率較低.Artman等(2010)提出逆時(shí)定位震源的自相關(guān)成像條件，得到震源能量團(tuán)，并以能量團(tuán)最大值處作為震源位置，但該方法在震源和檢波器之間存在大量干擾，影響定位分辨率；Nakata和Beroza(2016)對(duì)自相關(guān)成像算子進(jìn)行改進(jìn)，提出了互相關(guān)成像條件，互相關(guān)成像條件對(duì)地震記錄中的噪聲具有較好的壓制能力，可以得到具有較高空間分辨率的震源成像點(diǎn)，因而得到廣泛應(yīng)用；葛齊鑫等(2017)提出基于隨機(jī)分組的互相關(guān)成像條件，通過(guò)隨機(jī)分組進(jìn)行定位后再篩選成像結(jié)果來(lái)定位震源，在犧牲一部分計(jì)算效率的同時(shí)進(jìn)一步提高定位的抗噪性.Yuan等(2020)對(duì)最大振幅成像條件、自相關(guān)和互相關(guān)成像條件進(jìn)行多方面測(cè)試，比較了三種成像算法的運(yùn)行成本和成像分辨率，認(rèn)為在不同信噪比數(shù)據(jù)下應(yīng)平衡分辨率與抗噪性的需求，進(jìn)而優(yōu)選合適的成像算子.

在完全彈性介質(zhì)中，彈性波動(dòng)方程及聲波方程正演模擬對(duì)于時(shí)間反傳模擬是時(shí)不變的，通過(guò)完美的接收器采樣，反向傳播的波場(chǎng)會(huì)在時(shí)間上鏡像正向傳播的波場(chǎng).但是當(dāng)考慮地球介質(zhì)的固有衰減時(shí)，黏彈性波動(dòng)方程和黏聲波動(dòng)方程在時(shí)間反轉(zhuǎn)過(guò)程中不再是時(shí)不變的(Fink，2006)，逆時(shí)反傳過(guò)程中，反傳波場(chǎng)會(huì)進(jìn)一步衰減，即使接收器采樣完美，反向傳播的波場(chǎng)在時(shí)間上不會(huì)與正向傳播的波場(chǎng)對(duì)稱，并且重新聚焦的地震波能量耗散而且相位失真.為了恢復(fù)黏彈/黏聲介質(zhì)中反傳波場(chǎng)的時(shí)間對(duì)稱性，需要在逆時(shí)反傳過(guò)程中對(duì)耗散算子進(jìn)行改造以補(bǔ)償逆時(shí)反傳過(guò)程中的衰減，使反向傳播波場(chǎng)能夠恢復(fù)原始振幅并在準(zhǔn)確位置成像.Treeby和Cox(2010)在指數(shù)衰減黏聲介質(zhì)層析成像中討論了逆時(shí)反傳過(guò)程中的吸收補(bǔ)償問(wèn)題，對(duì)補(bǔ)償項(xiàng)進(jìn)行低通濾波以保持參數(shù)穩(wěn)定，改善了層析成像的總體分辨率；Zhu(2014)將解耦為頻散項(xiàng)和耗散項(xiàng)的黏聲波方程應(yīng)用到震源定位中，通過(guò)反轉(zhuǎn)耗散項(xiàng)達(dá)到補(bǔ)償衰減的作用，同時(shí)利用濾波消減補(bǔ)償引起的高頻噪聲問(wèn)題，實(shí)現(xiàn)了黏聲介質(zhì)下衰減補(bǔ)償?shù)恼鹪炊ㄎ唬玫搅溯^為準(zhǔn)確的震源空間位置；Zhu(2019)分離被動(dòng)源的波場(chǎng)并進(jìn)行了波場(chǎng)反傳，在不需要震源信息的情況下在層狀彈性介質(zhì)中實(shí)現(xiàn)了裂縫定位.Bai等(2019)對(duì)來(lái)自具有垂直對(duì)稱軸(VTI)的二維橫向各向同性介質(zhì)的合成微地震數(shù)據(jù)測(cè)試了衰減補(bǔ)償?shù)臅r(shí)間反轉(zhuǎn)成像算法，通過(guò)對(duì)多種因素與模型的分析驗(yàn)證了該方法適用于具有各向異性衰減的介質(zhì).Tang等(2020)對(duì)黏彈性正交各向異性介質(zhì)震源進(jìn)行各向異性衰減補(bǔ)償逆時(shí)定位，有效提高了水力壓裂監(jiān)測(cè)中震源空間位置定位精度.

本文基于Kjartansson常Q模型及Carcione的黏彈性波動(dòng)方程，通過(guò)low rank近似混合域算子提高計(jì)算效率，推導(dǎo)衰減補(bǔ)償?shù)酿椥越橘|(zhì)波動(dòng)方程及優(yōu)化成像算子，提高震源逆時(shí)定位的精度和計(jì)算效率.在含裂縫介質(zhì)中，利用模擬記錄與真實(shí)記錄的差異性分離散射波，基于衰減補(bǔ)償和分組互相關(guān)成像算子對(duì)裂縫進(jìn)行成像.

1 黏彈性介質(zhì)衰減補(bǔ)償

1.1 黏彈性介質(zhì)波動(dòng)方程

常Q模型(Kjartansson，1979)描述了一種品質(zhì)因子Q不隨頻率變化(但可以在空間上變化)的衰減介質(zhì)，衰減系數(shù)在頻率上是線性的.由常Q模型得到的松弛函數(shù)為：

(1)

其中M0是體積模量，Γ是伽馬函數(shù)，t0=1/ω0為參考時(shí)間，ω0為高于震源頻率的參考頻率，H是階躍函數(shù)，γ是介于0到1/2之間的無(wú)量綱參數(shù).

黏彈性各向同性介質(zhì)中的應(yīng)力-應(yīng)變(σ-ε)關(guān)系可以表示為：

(2)

式中符號(hào)*表示卷積，?2γε(t)/?t2γ表示應(yīng)變的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù).

當(dāng)Q→∞時(shí)，式(2)有完全彈性介質(zhì)本構(gòu)方程形式σ(t)=M0ε(t).對(duì)其傅里葉變換得到：

σ(ω)=M(ω)ε(ω),

(3)

式中ω表示角頻率，M(ω)=M0(iω/ω0)2γ表示復(fù)模量，根據(jù)復(fù)模量M與實(shí)模量M0的關(guān)系，相速度cP,S與衰減因子α表示為：

(4)

Carcione(2009)推導(dǎo)得出時(shí)間分?jǐn)?shù)階的黏彈性應(yīng)力-應(yīng)變關(guān)系：

(5)

圖1 (a) 速度與頻率關(guān)系； (b) 衰減系數(shù)與頻率關(guān)系Fig.1 (a) Relationship between velocity and frequency; (b) Relationship between attenuation coefficient and frequency

其中，δij是一個(gè)克羅內(nèi)克函數(shù)，i，j，k是空間指數(shù)，γ是介于0到1/2之間的無(wú)量綱參數(shù)，cP0,S0是在參考頻率ω0下的縱橫波速度，ρ為密度.

在二維情況下，式(5)可以表示為：

σxx=CλD2γP(εxx+εzz)-2CμD2γSεzz,

σzz=CλD2γP(εxx+εzz)-2CμD2γSεxx,

σxz=2CμD2γSεxz,

(6)

式(6)中的時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)可以近似為分?jǐn)?shù)階拉普拉斯算子通過(guò)傅里葉變換在波數(shù)域求出：

(-?2)γε=F-1{(k2)γF[ε]},

(-?2)γ-1/2ε=F-1{(k2)γ-1/2F[ε]},

(7)

其中F表示傅里葉變換，F(xiàn)-1表示反傅里葉變換.

綜上，在二維介質(zhì)條件下可以得到黏彈性波動(dòng)方程：

(8)

特別地，當(dāng)Q趨于∞時(shí)，γ趨于0，可得到彈性波動(dòng)方程.該式通過(guò)將時(shí)間分?jǐn)?shù)階導(dǎo)數(shù)轉(zhuǎn)移到空間波數(shù)域求解，僅需保存前一時(shí)刻波場(chǎng)，節(jié)省了內(nèi)存的同時(shí)提高計(jì)算效率，并且將介質(zhì)的頻散效應(yīng)與耗散效應(yīng)解耦，解耦項(xiàng)能夠在逆時(shí)延拓中進(jìn)行衰減補(bǔ)償.

1.2 Low rank近似

(9)

(9)式中空間域采用傅里葉計(jì)算精度較高，但是在時(shí)間域仍通過(guò)二階有限差分法計(jì)算，不可避免帶來(lái)一部分誤差，F(xiàn)irouzi等(2012)將sinc(vnΔtk/2)引入ikx得到新的波數(shù)-空間域混合矩陣，補(bǔ)償因時(shí)間差分帶來(lái)的誤差，得到：

(10)

設(shè)

W(x,k)=kxsinc(vPΔtk/2),

(11)

其中k表示波數(shù)向量，x表示空間向量.參考Fomel等(2013)的方法分解W(x，k)矩陣，簡(jiǎn)略步驟如下：

(1)從W矩陣中隨機(jī)選取n×R列矩陣組成WL0，n為期望得到的子矩陣秩數(shù)，R為采樣倍數(shù)，一般取4左右.對(duì)WL0進(jìn)行列主元正交三角分解，取其前n列組成候選矩陣WL1；

(2)與(1)同理隨機(jī)選取m×R列矩陣組成WR0，進(jìn)行行主元正交三角分解，取前m排組成候選矩陣WR1；

(3)定義一個(gè)系數(shù)矩陣an m，令W=WL1*an m*WR1，可通過(guò)矩陣運(yùn)算得到an m；

(4)合并三個(gè)子矩陣WL1*an m*WR1=W1，計(jì)算W1與W的誤差水平，若誤差在可接受范圍內(nèi)即得到W矩陣低秩近似形式，若誤差較高則從步驟1再次采樣計(jì)算，直到子矩陣符合誤差標(biāo)準(zhǔn).

通過(guò)上述分解，得到：

WM×N≈WL,M×nAn×mWR,m×N

(12)

其中n和m代表分解的秩數(shù)，秩數(shù)越高準(zhǔn)確性越高，但是相應(yīng)會(huì)增加計(jì)算量，n和m一般遠(yuǎn)小于M和N.WL由W(x，k)中的n列組成，與空間相關(guān)，WR由W(x，k)中的m行組成，與波數(shù)相關(guān).由此得到式(12)的低秩計(jì)算形式：

×F-1[WR(xm,k)ieikxΔx/2F(σxx)].

(13)

圖2 (a) 常規(guī)方法與low rank方法計(jì)算時(shí)間比較； (b) 解析解與數(shù)值解Fig.2 (a) Calculation time comparison between normal method and low rank method; (b) Analytical and numerical solutions

1.3 衰減補(bǔ)償逆時(shí)延拓

逆時(shí)反傳可以看作正演的逆過(guò)程，采用-t替換t，以黏聲介質(zhì)分?jǐn)?shù)階常Q波動(dòng)方程為例：

(14)

其中pB(xs,zs,t)=pF(xr,zr,T-t)，pB是時(shí)間t時(shí)的逆時(shí)空間波場(chǎng)，pF(xr,zr,T-t)是接收器位置在T-t時(shí)刻的數(shù)據(jù)記錄.方程右側(cè)兩項(xiàng)分別為解耦的頻散項(xiàng)和耗散項(xiàng).頻散項(xiàng)與時(shí)間無(wú)關(guān)，而耗散項(xiàng)隨時(shí)間會(huì)使得地震波能量逐漸減弱，故在地震波反傳過(guò)程中應(yīng)該對(duì)波場(chǎng)進(jìn)行衰減補(bǔ)償(Zhu and Carcione，2014)，即將耗散項(xiàng)顛倒符號(hào)得到式(15)：

(15)

由正演計(jì)算得到的衰減與頻率成線性正比關(guān)系，在波場(chǎng)正傳過(guò)程中，能量逐漸耗散可以保持穩(wěn)定，但是反向延拓過(guò)程中由于補(bǔ)償項(xiàng)的存在，高頻能量得到極大增強(qiáng)，而地震記錄中高頻信號(hào)可能會(huì)被噪聲污染，補(bǔ)償過(guò)程會(huì)將這些噪聲增強(qiáng)導(dǎo)致波動(dòng)方程不穩(wěn)定(Treeby et al.，2010)，所以需要對(duì)逆時(shí)延拓中的波場(chǎng)進(jìn)行低通濾波，截止波數(shù)可以通過(guò)介質(zhì)最大速度的截?cái)囝l率計(jì)算獲得.截?cái)囝l率需根據(jù)地震記錄中有效信號(hào)和噪聲的頻率分布確定：圖3a為截?cái)囝l率為100 Hz情況下四階巴特沃茲低通濾波器示例，在截?cái)囝l率處濾波系數(shù)為0.707，該濾波器有平穩(wěn)的幅頻特性；圖3b為選取試驗(yàn)過(guò)程中記錄數(shù)據(jù)的單個(gè)檢波器數(shù)據(jù)進(jìn)行分析得到的頻譜(紅色線)及多個(gè)檢波器計(jì)算的平均振幅譜(黑色線)，可以看到數(shù)據(jù)段有效信息大多集中在70 Hz以下，為保護(hù)波場(chǎng)的有效信息，將截?cái)囝l率稍增大，設(shè)置在100 Hz限制補(bǔ)償?shù)念l率段來(lái)穩(wěn)定波場(chǎng)，圖3c為對(duì)應(yīng)波數(shù)域低通濾波器示例.

圖3 (a) 濾波器響應(yīng)，虛線處為截?cái)囝l率對(duì)應(yīng)的濾波系數(shù)； (b) 含噪數(shù)據(jù)段頻譜及整體平均振幅譜分析； (c) 波數(shù)域?yàn)V波器Fig.3 (a) Filter response. The dotted line is the filter coefficient of cutoff frequency; (b) Analysis of the frequency spectrum and overall average amplitude spectrum of a noisy data; (c) Wavenumber domain filter

在對(duì)波場(chǎng)進(jìn)行濾波時(shí)，不應(yīng)對(duì)傳播算子信號(hào)進(jìn)行處理，需要將傳播算子從頻散項(xiàng)中分離(Zhu，2015)，對(duì)式(15)進(jìn)行改造得到：

(16)

根據(jù)式(16)的思想，對(duì)式(8)進(jìn)行改造得到黏彈性介質(zhì)逆時(shí)反傳波動(dòng)方程：

σxx=[Cλ(εxx+εzz)-2Cμεzz]+FLP[ηP(-?2)γP(εxx

+εzz)-2ηS(-?2)γSεzz-Cλ(εxx+εzz)+2Cμεzz]

σzz=[Cλ(εxx+εzz)-2Cμεxx]+FLP[ηP(-?2)γP(εxx

+εzz)-2ηS(-?2)γSεxx-Cλ(εxx+εzz)+2Cμεxx]

σxz=2Cμεxz+FLP[2ηS(-?2)γSεxz-2Cμεxz]

(17)

FLP代表低通濾波器，本文選用巴特沃茲低通濾波器，在波數(shù)域其表達(dá)式為：

(18)

其中f表示頻率，fc表示截?cái)囝l率，n代表濾波器階數(shù).

在黏彈介質(zhì)中實(shí)現(xiàn)逆時(shí)反傳的過(guò)程包括三個(gè)步驟：

(1)將地震記錄逆時(shí)顛倒，然后將該數(shù)據(jù)作為原始接收器位置的邊界條件；

(2)使邊界波在相反的時(shí)間內(nèi)傳播并濾波，隨著時(shí)間的流逝，耗散和頻散效應(yīng)會(huì)自動(dòng)得到補(bǔ)償；

(3)通過(guò)適當(dāng)?shù)某上駰l件搜索震源位置.

2 微地震衰減補(bǔ)償逆時(shí)定位

2.1 衰減補(bǔ)償效果分析

設(shè)置二維介質(zhì)中震源、接收點(diǎn)分布如圖4a，對(duì)2號(hào)環(huán)狀檢波器組接收到的地震記錄進(jìn)行逆時(shí)延拓，得到近震源1號(hào)檢波器處波場(chǎng)如圖4b，在0.12 s之前的噪聲是反傳中互相串?dāng)_的信號(hào)以及前文提到的隨機(jī)噪聲，可通過(guò)成像條件進(jìn)行壓制.放大0.12～0.24 s部分如圖4c，可以看到真實(shí)波場(chǎng)(實(shí)線VE)與反傳波場(chǎng)(VE-VE)存在一定能量損失，其原因是接收點(diǎn)不能捕獲所有能量，而且為了壓制噪聲需要對(duì)信號(hào)進(jìn)行一定程度低通濾波，也會(huì)造成少量能量損失.圖中未補(bǔ)償波場(chǎng)(VE-EL)明顯出現(xiàn)延滯，而進(jìn)行衰減補(bǔ)償后波場(chǎng)傳播時(shí)間已經(jīng)校正.通過(guò)衰減補(bǔ)償?shù)膌ow rank波動(dòng)方程進(jìn)行波場(chǎng)延拓可以在正確位置成像，修正黏彈性介質(zhì)對(duì)地震波的吸收衰減作用.

圖4 (a) 震源與環(huán)形檢波器排列； (b) 逆時(shí)延拓波場(chǎng)對(duì)比； (c) 0.12～0.24 s波場(chǎng)對(duì)比Fig.4 (a) Distribution of source and circle receiver; (b) Comparison of backpropagation wavefield; (c) Wavefield comparison at 0.12～0.24 s

2.2 逆時(shí)定位效果比較

設(shè)置一個(gè)簡(jiǎn)單的三層模型，大小為2000 m×2000 m，網(wǎng)格間距10 m，參數(shù)見(jiàn)表1，在(1000 m，1300 m)處正應(yīng)力分量加載主頻30 Hz雷克子波，分別在地表和700 m處一1500 m深井記錄地震數(shù)據(jù).得到多分量地震記錄后(如圖5所示)，可以進(jìn)行單分量反傳定位或多分量反傳定位，地表地震記錄與井中地震記錄聯(lián)合反傳定位，也可進(jìn)行分頻段多尺度反傳(本文默認(rèn)進(jìn)行全頻段記錄反傳)等各種延拓方法.首先在黏彈性介質(zhì)中分別進(jìn)行彈性反傳(VE-EL，圖6a、d)、無(wú)補(bǔ)償衰減反傳(VE-QE，圖6b、e)和衰減補(bǔ)償反傳(VE-VE，圖6c、f)，并通過(guò)自相關(guān)成像算子進(jìn)行震源定位，地面定位結(jié)果提取地表1000 m處應(yīng)力垂直分量如圖6g，通過(guò)衰減補(bǔ)償定位震源深度為1300～1310 m，無(wú)補(bǔ)償時(shí)定位在1260～1280 m，使用彈性波動(dòng)方程定位在1250～1280 m處.無(wú)補(bǔ)償波能量分布在地表到震源之間，分辨率較低，彈性波能量聚集性弱于衰減補(bǔ)償波場(chǎng)，通過(guò)衰減補(bǔ)償后波場(chǎng)不但定位誤差最小，且能量聚集更集中便于識(shí)別；井中地震數(shù)據(jù)提取井中1300 m處數(shù)據(jù)應(yīng)力水平分量圖6h，震源定位效果與地面定位效果基本相同，但是在等深度處400 m外出現(xiàn)另一峰值，這是由于探測(cè)井的位置分布及井深的局限性，導(dǎo)致有效信息空間采樣不足，殘余能量聚集到井對(duì)稱位置形成假震源，由于其并非真實(shí)震源，故波場(chǎng)能量弱于真實(shí)震源，在復(fù)雜介質(zhì)中假震源聚焦能力會(huì)進(jìn)一步減弱，且補(bǔ)償后假震源能量更小有助于識(shí)別.通過(guò)染色算法(Chen and Jia,2014)或地面記錄與井中記錄聯(lián)合反傳(圖6i)等方法可以規(guī)避假震源的影響.

表1 三層模型參數(shù)Table 1 Three-layer model parameters

圖5 (a) 地表水平分量記錄； (b) 地表垂直分量記錄； (c) 井中水平分量記錄； (d) 井中垂直分量記錄Fig.5 (a) The horizontal component of the surface record; (b) The vertical component of the surface record; (c) The horizontal component of the borehole record; (d) The vertical component of the borehole record

圖6 地面接收情況下(a)彈性介質(zhì), (b)黏彈性介質(zhì)無(wú)補(bǔ)償和(c)黏彈性介質(zhì)衰減補(bǔ)償定位結(jié)果; 井中接收情況下(d)彈性介質(zhì), (e)黏彈性介質(zhì)無(wú)補(bǔ)償和(f)黏彈性介質(zhì)衰減補(bǔ)償定位結(jié)果; (g)地面接收定位結(jié)果單道數(shù)據(jù)對(duì)比，虛線為三種方法定位震源位置; (h)井中接收定位結(jié)果單道數(shù)據(jù)對(duì)比，虛線為三種方法定位震源位置; (i)地面與井中記錄衰減補(bǔ)償聯(lián)合定位，虛線交點(diǎn)為真實(shí)震源位置Fig.6 Surface monitoring source location result of (a) elastic media, (b) viscoelastic media without attenuation compensation and (c) viscoelastic media with attenuation compensation; Borehole monitoring source location result of (d) elastic media, (e) viscoelastic media without attenuation compensation and (f) viscoelastic media with attenuation compensation; (g) Surface monitoring source location comparison of single channel data. The dotted lines are source location results of three methods; (h) Borehole monitoring source location comparison of single channel data. The dotted lines are source location results of three methods; (i) Surface and borehole record joint location result with attenuation compensation. The intersection of dotted lines is the actual source location

2.3 成像算子優(yōu)化

圖7a為Marmousi模型，模型大小為3000 m×3000 m，網(wǎng)格間距10 m，將震源設(shè)置在(1500 m，1500 m)處，在地表采集地震信號(hào)，地表檢波器道間距100 m，在震源定位過(guò)程中需要考慮各種干擾，如采集信號(hào)中存在的噪聲、地下介質(zhì)參數(shù)的不確定性等等.將Marmousi模型中的縱波速度、密度參數(shù)進(jìn)行平滑(如圖7b、c)，設(shè)介質(zhì)泊松比為0.25，計(jì)算得到橫波速度，縱波品質(zhì)因子設(shè)為縱波速度的1/40，橫波品質(zhì)因子設(shè)為縱波品質(zhì)因子的0.83倍.在地震記錄中加入隨機(jī)噪聲，進(jìn)行自相關(guān)、互相關(guān)及結(jié)合兩種方法的優(yōu)化互相關(guān)算法進(jìn)行震源定位.地震記錄加入強(qiáng)噪聲后如圖8，含噪信號(hào)信噪比約為-18 dB，有效信號(hào)幾乎被噪聲淹沒(méi).

圖7 (a) Marmousi縱波速度模型； (b) 平滑后Marmousi縱波速度； (c) 平滑后的密度模型Fig.7 (a) Marmousi P-wave velocity model; (b) Smoothed Marmousi P-wave velocity model; (c) Smoothed density model

圖8 含噪的地表接收數(shù)據(jù)(a) 水平速度分量； (b) 垂直速度分量.Fig.8 Noisy surface data(a) Horizontal velocity component; (b) Vertical velocity component.

成像算子中，最大值成像法(McMechan，1982)是在逆時(shí)反傳過(guò)程中不斷檢索波場(chǎng)能量最大值，以最終的最大值位置為震源空間位置，由于多次檢索導(dǎo)致計(jì)算效率較低；自相關(guān)成像算子(Artman et al.，2010)在波場(chǎng)延拓過(guò)程中進(jìn)行零延時(shí)自相關(guān)，對(duì)背景噪聲有一定的壓制能力，但是對(duì)記錄中的隨機(jī)噪聲無(wú)法進(jìn)行有效壓制(圖9a)；互相關(guān)成像算子(Nakata and Beroza，2016)利用噪聲之間相干性較低的特點(diǎn)，將各檢波點(diǎn)記錄單獨(dú)進(jìn)行逆時(shí)延拓并進(jìn)行互相關(guān)，壓制噪聲能力較強(qiáng)，而隨檢波器的增加計(jì)算成本呈線性增長(zhǎng)，計(jì)算效率和運(yùn)行內(nèi)存要額外消耗數(shù)十倍至上百倍之多，為節(jié)約計(jì)算成本一般按位置分組逆時(shí)延拓.圖9b為將檢波器分為等段距離的三組后進(jìn)行互相關(guān)得到的定位圖像，雖然壓制噪聲能力較自相關(guān)算法更強(qiáng)，但是其能量最大值位于模型上方一斷層界面處.圖9c為在等距三組檢波器分組的基礎(chǔ)上針對(duì)成像算子進(jìn)行優(yōu)化，在互相關(guān)算法的基礎(chǔ)上進(jìn)行一次自相關(guān)，結(jié)合自相關(guān)成像算子對(duì)背景噪聲的壓制和互相關(guān)成像算子對(duì)隨機(jī)噪聲壓制的優(yōu)點(diǎn)，加強(qiáng)成像算子的聚焦能力(Tang et al.，2020).

圖9 (a) 自相關(guān)成像； (b) 三組檢波器互相關(guān)成像； (c) 三組檢波器優(yōu)化互相關(guān)成像，虛線交點(diǎn)為真實(shí)震源位置Fig.9 (a) Autocorrelation imaging; (b) Cross-correlation imaging with three groups of receivers; (c) Optimized cross-correlation imaging with three groups of receivers. The intersection of dotted lines is actual source location

對(duì)于互相關(guān)算子中檢波器的分組，傳統(tǒng)檢波器陣列為[1,1…1,2,2…2,…n,n…n]，將檢波器按位置分組，壓制n組之間的不相干噪聲，當(dāng)n較小時(shí)成像效果可能不佳，應(yīng)對(duì)方法是提高分組數(shù)量，這種做法雖然會(huì)提高成像效果，但帶來(lái)的計(jì)算成本違背了分組的初衷.葛奇鑫等(2017)提出隨機(jī)分組方法，由于需要多次分組成像再進(jìn)行篩選震源，壓制噪聲較好但是效率相比前者要低，本文不做討論對(duì)比.通過(guò)將檢波器陣列穿插劃分成[1,2,3,…,n,1,2,3,…n…1,2,3…n]，在不提高計(jì)算量的情況下得到逆時(shí)成像如圖10a，在減少一個(gè)分組時(shí)兩種分組方法得到成像如圖10b、c .

圖10 (a) 三組檢波器穿插陣列成像； (b) 兩組檢波器常規(guī)互相關(guān)成像； (c) 兩組檢波器穿插陣列成像，其中虛線交點(diǎn)為真實(shí)震源位置Fig.10 (a) Interleaved array imaging with three groups of receivers; (b) Conventional cross-correlation imaging with two groups of receivers; (c) Interleaved array imaging with two groups of receivers. The intersection of dotted lines is the actual source location

對(duì)比圖10a和圖9c，同樣計(jì)算量下，穿插陣列分組能夠更好地壓制噪聲，而在減少1/3計(jì)算量情況下，常規(guī)分組成像定位結(jié)果向左方偏離，而穿插陣列分組能夠得到與原三個(gè)分組接近的定位精度.穿插陣列互相關(guān)算子能夠在噪聲壓制效果和降低計(jì)算成本上獲得更好的效果.改變模型參數(shù)平滑尺度(圖11a、b)，進(jìn)行定位的結(jié)果如圖11c、d，模型參數(shù)的平滑會(huì)導(dǎo)致能量聚焦性減弱.在速度整體變化時(shí)，由于本方法基于波場(chǎng)能量進(jìn)行定位，無(wú)法克服速度變化帶來(lái)的干擾，可以結(jié)合地表地震數(shù)據(jù)通過(guò)全波形反演以獲得更好的速度模型.

圖11 (a)和(b)為不同平滑程度的速度模型； (c) a模型定位結(jié)果； (d) b模型定位結(jié)果，其中虛線交點(diǎn)為真實(shí)震源位置Fig.11 (a) and (b) are velocity models with different level of smoothness; (c) The location result of a model; (d) The location result of b model. The intersection of dotted lines is the actual source location

2.4 三維模型測(cè)試

通過(guò)三維overthrust模型進(jìn)行定位測(cè)試，模型大小1500 m×1500 m×1200 m，模型縱波速度參數(shù)如圖12a所示，在地表設(shè)置四條測(cè)線進(jìn)行接收(圖12b)，爆炸源設(shè)置在測(cè)線交點(diǎn)下方深960 m處.圖12d為地表接收的592道垂直分量記錄，在地表記錄中加入信噪比為1.25 dB噪聲，并通過(guò)平滑后的模型(圖12c)進(jìn)行定位.定位結(jié)果如圖13，其中圖13a為通過(guò)衰減補(bǔ)償定位的結(jié)果，圖13b是未補(bǔ)償定位結(jié)果，圖13c為在真實(shí)震源處抽取的垂直單道歸一化數(shù)據(jù)對(duì)比，虛線位置為真實(shí)震源深度.從圖13a、b可以看到衰減補(bǔ)償能量聚焦在震源附近，而未補(bǔ)償定位結(jié)果聚焦效果較差，無(wú)法辨識(shí)震源；圖13c中補(bǔ)償后震源定位結(jié)果準(zhǔn)確，而未進(jìn)行補(bǔ)償?shù)亩ㄎ唤Y(jié)果其能量集中在地表低速帶附近，震源處能量峰值位置也偏淺，存在一定誤差.如果提高平滑度，即模型參數(shù)不準(zhǔn)確時(shí)，通過(guò)衰減補(bǔ)償也可能無(wú)法定位震源真實(shí)位置，因此獲得準(zhǔn)確的初始速度模型對(duì)震源定位的效果影響較高.本文方法也可以用于各向異性黏彈性介質(zhì)，對(duì)于各向異性黏彈性模型的衰減補(bǔ)償逆時(shí)定位結(jié)果見(jiàn)附錄A.

圖12 (a) Overthrust模型； (b) 檢波器及震源位置分布； (c) 用于定位的平滑速度模型； (d) 微地震記錄Fig.12 (a) Overthrust model; (b) Distribution of receivers and source; (c) The smoothed velocity model for location; (d) Microseismic record

圖13 (a) 衰減補(bǔ)償定位結(jié)果； (b) 無(wú)補(bǔ)償定位； (c) 垂直單道數(shù)據(jù)對(duì)比，虛線對(duì)應(yīng)真實(shí)震源位置Fig.13 (a) Location result with attenuation compensation; (b) Location result without attenuation compensation; (c) Comparison of vertical single channel data. The dotted line corresponds to the actual source location

3 裂縫成像

在水力壓裂過(guò)程中，地下存在裂縫時(shí)，根據(jù)惠更斯原理，地震波經(jīng)過(guò)裂縫時(shí)會(huì)產(chǎn)生透射波和散射波，散射波可以視為以散射點(diǎn)(即裂縫)為震源發(fā)出的子波，如果能將散射波與透射波分離，對(duì)散射波和透射波進(jìn)行互相關(guān)反傳，它們將會(huì)在裂縫處聚焦，實(shí)現(xiàn)被動(dòng)源裂縫成像.

3.1 井中接收數(shù)據(jù)裂縫成像

設(shè)置一個(gè)層狀模型對(duì)裂縫成像進(jìn)行模擬，縱波速度模型如圖14a，模型分為三層，大小為600 m×600 m，網(wǎng)格間距2 m，在第二層存在兩條裂縫，模型參數(shù)如表2，在(100 m，2800 m)處設(shè)置微地震震源，并于地表接收地震記錄，其中水平速度分量如圖15a.對(duì)圖14a中的模型參數(shù)進(jìn)行平滑得到的縱波速度模型如圖14b所示.

表2 裂縫模型參數(shù)Table 2 Fracture model parameters

圖14 (a) 裂縫模型，其中“*”為震源，為檢波器； (b) 平滑模型Fig.14 (a) Fracture model. ‘*’ is the source, and is the receiver; (b) Smoothed model

通過(guò)檢索每一道地震記錄的走時(shí)，將初至波拉平處理；然后對(duì)地震記錄進(jìn)行頻譜分析，通過(guò)主頻與時(shí)間采樣間隔估算地震波波長(zhǎng)；再按照波長(zhǎng)切割分離初至透射波和散射波，并通過(guò)前文介紹衰減補(bǔ)償逆時(shí)反傳方法進(jìn)行裂縫成像，得到成像結(jié)果如圖15d，井中接收數(shù)據(jù)逆時(shí)波場(chǎng)能量聚集在裂縫的邊緣位置，誤差較小.

圖15 (a) 井中vx分量記錄； (b) 透射波記錄； (c) 散射波及層間反射波記錄； (d) 裂縫成像結(jié)果Fig.15 (a) Borehole record of vx component; (b) Record of transmitted wave; (c) Records of scattering wave and reflected wave; (d) Fracture imaging result

3.2 地面接收數(shù)據(jù)裂縫成像

在地表進(jìn)行裂縫定位時(shí)，受層間多次波的影響，不能簡(jiǎn)單地切除初至波進(jìn)行散射波分離，本文通過(guò)模擬新的波場(chǎng)并進(jìn)行處理達(dá)到分離散射波的目的.設(shè)置縱波速度模型如圖16a，模型分為四層，大小為500 m×250 m，網(wǎng)格間距1 m，在第三層存在四條裂縫，模型參數(shù)如表3，在(250 m，240 m)處設(shè)置微地震震源，并于地表接收地震記錄，其中垂直速度分量如圖17a.對(duì)圖16a中的速度模型進(jìn)行平滑得到的速度模型如圖16b所示.

表3 裂縫模型參數(shù)Table 3 Fracture model parameters

圖16 (a) 裂縫模型，其中“*”為震源，為檢波器； (b) 平滑模型； (c) 震源定位結(jié)果，其中虛線交點(diǎn)為真實(shí)震源位置，“*”為定位結(jié)果位置Fig.16 (a) Fracture model; “*” is the source, and is the receiver; (b) Smoothed model; (c) Source location result. The intersection of dotted lines is the actual source location, and “*” is source location result

(1)首先通過(guò)衰減補(bǔ)償逆時(shí)定位，尋找微地震震源的空間位置，圖16c中“*”為定位結(jié)果，虛線交點(diǎn)為真實(shí)震源位置，可以看到定位結(jié)果準(zhǔn)確.

(2)在定位震源位置處激發(fā)震源，對(duì)圖16b模型進(jìn)行二次波場(chǎng)模擬，得到模擬記錄垂直速度分離如圖17b.由于已知參數(shù)不夠精準(zhǔn)，模擬直達(dá)波與真實(shí)直達(dá)波會(huì)存在走時(shí)和振幅的偏差.

(3)記錄檢波器每一道的最大值，將所有記錄按道進(jìn)行歸一化處理.通過(guò)檢索真實(shí)記錄直達(dá)波走時(shí)，將模擬記錄的直達(dá)波校正至真實(shí)直達(dá)波位置處，然后將兩個(gè)記錄進(jìn)行相減，并將新得到的記錄進(jìn)行逆歸一化，每道進(jìn)行同樣處理后得到散射波如圖17c.

(4)對(duì)分離后的散射波記錄與處理后的模擬波記錄進(jìn)行互相關(guān)反傳，并切除檢測(cè)到的震源能量團(tuán)后成像如圖17e，其中右三處裂縫位置成像能量較強(qiáng)，左側(cè)裂縫由于距震源最遠(yuǎn)，成像能量較弱.

另外應(yīng)用Zhu(2019)提出的中值濾波方法進(jìn)行測(cè)試，得到散射波如圖17d，成像結(jié)果如圖17f，其結(jié)果與本文方法結(jié)果相似，左側(cè)裂縫能量較弱.

圖17 (a) 裂縫模型合成記錄vz分量； (b) 平滑模型合成記錄vz分量； (c) 匹配分離散射波； (d) 中值濾波分離散射波； (e) 匹配分離散射波裂縫成像； (f) 中值濾波分離散射波裂縫成像Fig.17 (a) Synthetic record vz component of fracture model; (b) Synthetic record vz component of smoothed model; (c) Match-separated scattering wave; (d) Median filter separated scattering wave; (e) The fracture imaging of match-separated scattering wave; (f) The fracture imaging of median filter separated scattering wave

4 結(jié)論

針對(duì)黏彈性介質(zhì)震源定位及裂縫成像中計(jì)算效率及準(zhǔn)確度的問(wèn)題，本文基于Kjartansson常Q模型及Carcione的黏彈性介質(zhì)波動(dòng)方程，引入low rank算法進(jìn)行計(jì)算加速，并對(duì)成像算子進(jìn)行了改進(jìn)，實(shí)現(xiàn)衰減補(bǔ)償震源定位及裂縫成像，得到以下結(jié)論和認(rèn)識(shí)：

(1)通過(guò)分?jǐn)?shù)階拉普拉斯算子將能量耗散和相位頻散明確分離，在逆時(shí)反傳時(shí)通過(guò)反轉(zhuǎn)耗散算子的符號(hào)并使保持頻散項(xiàng)符號(hào)不變來(lái)完成衰減補(bǔ)償.補(bǔ)償后的波場(chǎng)能夠恢復(fù)能量并校正能量聚集位置，定位震源真實(shí)位置.

(2)針對(duì)計(jì)算混合域算子時(shí)逐點(diǎn)傅里葉變換計(jì)算成本較高的問(wèn)題，通過(guò)low rank分解將混合域算子分解為三組矩陣，采用巴特沃茲低通濾波器穩(wěn)定衰減補(bǔ)償波場(chǎng)，避免了高頻噪聲被補(bǔ)償引起的高能量波場(chǎng)掩蓋有效信號(hào)，在保證波場(chǎng)模擬精度的情況下降低了計(jì)算成本.

(3)針對(duì)逆時(shí)定位過(guò)程中存在的隨機(jī)噪聲影響，優(yōu)化成像算子，在減少計(jì)算量的同時(shí)提高定位分辨率.優(yōu)化后的成像算子能夠在加入強(qiáng)隨機(jī)噪聲的地震記錄中高效率高精度定位震源位置.

(4)對(duì)于水力壓裂過(guò)程中介質(zhì)中存在的裂縫，通過(guò)分離地震記錄中的透射波和散射波并進(jìn)行互相關(guān)反傳，能夠?qū)α芽p進(jìn)行成像，但實(shí)際上常規(guī)方法對(duì)散射波的分離效果仍不理想，新近發(fā)展的深度學(xué)習(xí)方法可能是一個(gè)解決方向.

由于采用常規(guī)的波場(chǎng)能量作為成像值，在速度存在整體偏離的情況下震源位置存在誤差.對(duì)于速度整體的偏離引起的誤差可以與全波形反演方法結(jié)合，本文采用的成像算子也可以移植到能量范數(shù)、反褶積算法等其他算法上.震源定位問(wèn)題涉及初始模型參數(shù)、各向異性、黏彈性以及震源機(jī)制等多方面，本文主要解決黏彈性介質(zhì)帶來(lái)的地震波衰減的影響，并未對(duì)其他因素詳細(xì)討論，涉及多參數(shù)影響下的震源定位問(wèn)題，或許可以借助全波形反演的思想進(jìn)行解決，尚需進(jìn)一步研究.

附錄A 三維黏彈性各向異性介質(zhì)震源定位

本文討論了二維各向同性黏彈性介質(zhì)的衰減補(bǔ)償逆時(shí)定位，該方法在各向異性介質(zhì)中同樣適用.這里設(shè)計(jì)了一個(gè)黏彈性VTI介質(zhì)進(jìn)行測(cè)試，對(duì)離散采樣的overthrust模型(如圖A1a)添加速度各向異性，設(shè)置Thomsen參數(shù)ε=0.2，δ=0.1，γ=0.15，震源為加載在中心點(diǎn)深900 m處的爆炸源，地表設(shè)置四條測(cè)線，利用平滑后的模型進(jìn)行逆時(shí)定位，成像算子采用本文設(shè)計(jì)的優(yōu)化互相關(guān)算子.圖A1b為檢波器及震源位置，圖A1c為地表592道垂直分量記錄，圖A1d為震源定位結(jié)果剖面，圖A1e為抽取震源位置垂直單道歸一化曲線，可以看到，在模型中添加速度各向異性并不影響本文方法的定位效果，在模型參數(shù)較準(zhǔn)確時(shí)，本文方法能夠準(zhǔn)確定位震源.

圖A1 (a) Overthrust模型； (b) 檢波器及震源位置； (c) 地表記錄； (d) 震源定位結(jié)果剖面；(e) 定位結(jié)果垂直單道數(shù)據(jù)，虛線對(duì)應(yīng)真實(shí)震源深度Fig.A1 (a) Overthrust model; (b) Receivers and source location; (c) Surface record; (d) Source location result profile; (e) Vertical single-channel data of source location. The dotted line corresponds to the actual source depth