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

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

0 引言

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

(3)

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

(4)

Carcione(2009)推導(dǎo)得出時間分?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是一個克羅內(nèi)克函數(shù)，i，j，k是空間指數(shù)，γ是介于0到1/2之間的無量綱參數(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)中的時間分?jǐn)?shù)階導(dǎo)數(shù)可以近似為分?jǐn)?shù)階拉普拉斯算子通過傅里葉變換在波數(shù)域求出：

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

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

(7)

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

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

(8)

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

1.2 Low rank近似

(9)

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

(10)

設(shè)

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

(11)

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

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

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

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

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

通過上述分解，得到：

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

(12)

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

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

(13)

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

1.3 衰減補償逆時延拓

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

(14)

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

(15)

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

圖3 (a) 濾波器響應(yīng)，虛線處為截斷頻率對應(yīng)的濾波系數(shù)； (b) 含噪數(shù)據(jù)段頻譜及整體平均振幅譜分析； (c) 波數(shù)域濾波器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

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

(16)

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

σ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表示截斷頻率，n代表濾波器階數(shù).

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

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

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

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

2 微地震衰減補償逆時定位

2.1 衰減補償效果分析

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

圖4 (a) 震源與環(huán)形檢波器排列； (b) 逆時延拓波場對比； (c) 0.12～0.24 s波場對比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è)置一個簡單的三層模型，大小為2000 m×2000 m，網(wǎng)格間距10 m，參數(shù)見表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)、無補償衰減反傳(VE-QE，圖6b、e)和衰減補償反傳(VE-VE，圖6c、f)，并通過自相關(guān)成像算子進(jìn)行震源定位，地面定位結(jié)果提取地表1000 m處應(yīng)力垂直分量如圖6g，通過衰減補償定位震源深度為1300～1310 m，無補償時定位在1260～1280 m，使用彈性波動方程定位在1250～1280 m處.無補償波能量分布在地表到震源之間，分辨率較低，彈性波能量聚集性弱于衰減補償波場，通過衰減補償后波場不但定位誤差最小，且能量聚集更集中便于識別；井中地震數(shù)據(jù)提取井中1300 m處數(shù)據(jù)應(yīng)力水平分量圖6h，震源定位效果與地面定位效果基本相同，但是在等深度處400 m外出現(xiàn)另一峰值，這是由于探測井的位置分布及井深的局限性，導(dǎo)致有效信息空間采樣不足，殘余能量聚集到井對稱位置形成假震源，由于其并非真實震源，故波場能量弱于真實震源，在復(fù)雜介質(zhì)中假震源聚焦能力會進(jìn)一步減弱，且補償后假震源能量更小有助于識別.通過染色算法(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ì)無補償和(c)黏彈性介質(zhì)衰減補償定位結(jié)果; 井中接收情況下(d)彈性介質(zhì), (e)黏彈性介質(zhì)無補償和(f)黏彈性介質(zhì)衰減補償定位結(jié)果; (g)地面接收定位結(jié)果單道數(shù)據(jù)對比，虛線為三種方法定位震源位置; (h)井中接收定位結(jié)果單道數(shù)據(jù)對比，虛線為三種方法定位震源位置; (i)地面與井中記錄衰減補償聯(lián)合定位，虛線交點為真實震源位置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)處，在地表采集地震信號，地表檢波器道間距100 m，在震源定位過程中需要考慮各種干擾，如采集信號中存在的噪聲、地下介質(zhì)參數(shù)的不確定性等等.將Marmousi模型中的縱波速度、密度參數(shù)進(jìn)行平滑(如圖7b、c)，設(shè)介質(zhì)泊松比為0.25，計算得到橫波速度，縱波品質(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，含噪信號信噪比約為-18 dB，有效信號幾乎被噪聲淹沒.

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

圖9 (a) 自相關(guān)成像； (b) 三組檢波器互相關(guān)成像； (c) 三組檢波器優(yōu)化互相關(guān)成像，虛線交點為真實震源位置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

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

圖10 (a) 三組檢波器穿插陣列成像； (b) 兩組檢波器常規(guī)互相關(guān)成像； (c) 兩組檢波器穿插陣列成像，其中虛線交點為真實震源位置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

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

圖11 (a)和(b)為不同平滑程度的速度模型； (c) a模型定位結(jié)果； (d) b模型定位結(jié)果，其中虛線交點為真實震源位置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 三維模型測試

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

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

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

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

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

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

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

圖16 (a) 裂縫模型，其中“*”為震源，為檢波器； (b) 平滑模型； (c) 震源定位結(jié)果，其中虛線交點為真實震源位置，“*”為定位結(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)首先通過衰減補償逆時定位，尋找微地震震源的空間位置，圖16c中“*”為定位結(jié)果，虛線交點為真實震源位置，可以看到定位結(jié)果準(zhǔn)確.

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

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

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

另外應(yīng)用Zhu(2019)提出的中值濾波方法進(jìn)行測試，得到散射波如圖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é)論

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

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

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

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

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

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

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

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

圖A1 (a) Overthrust模型； (b) 檢波器及震源位置； (c) 地表記錄； (d) 震源定位結(jié)果剖面；(e) 定位結(jié)果垂直單道數(shù)據(jù)，虛線對應(yīng)真實震源深度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