基于Low rank近似的黏彈性介質衰減補償微地震逆時定位與裂縫成像

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

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

0 引言

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

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

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

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

本文基于Kjartansson常Q模型及Carcione的黏彈性波動方程，通過low rank近似混合域算子提高計算效率，推導衰減補償的黏彈性介質波動方程及優化成像算子，提高震源逆時定位的精度和計算效率.在含裂縫介質中，利用模擬記錄與真實記錄的差異性分離散射波，基于衰減補償和分組互相關成像算子對裂縫進行成像.

1 黏彈性介質衰減補償

1.1 黏彈性介質波動方程

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

(1)

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

黏彈性各向同性介質中的應力-應變(σ-ε)關系可以表示為：

(2)

式中符號*表示卷積，?2γε(t)/?t2γ表示應變的時間分數階導數.

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

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

(3)

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

(4)

Carcione(2009)推導得出時間分數階的黏彈性應力-應變關系：

(5)

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

其中，δij是一個克羅內克函數，i，j，k是空間指數，γ是介于0到1/2之間的無量綱參數，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)中的時間分數階導數可以近似為分數階拉普拉斯算子通過傅里葉變換在波數域求出：

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

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

(7)

其中F表示傅里葉變換，F-1表示反傅里葉變換.

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

(8)

特別地，當Q趨于∞時，γ趨于0，可得到彈性波動方程.該式通過將時間分數階導數轉移到空間波數域求解，僅需保存前一時刻波場，節省了內存的同時提高計算效率，并且將介質的頻散效應與耗散效應解耦，解耦項能夠在逆時延拓中進行衰減補償.

1.2 Low rank近似

(9)

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

(10)

設

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

(11)

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

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

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

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

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

通過上述分解，得到：

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

(12)

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

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

(13)

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

1.3 衰減補償逆時延拓

逆時反傳可以看作正演的逆過程，采用-t替換t，以黏聲介質分數階常Q波動方程為例：

(14)

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

(15)

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

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

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

(16)

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

σ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代表低通濾波器，本文選用巴特沃茲低通濾波器，在波數域其表達式為：

(18)

其中f表示頻率，fc表示截斷頻率，n代表濾波器階數.

在黏彈介質中實現逆時反傳的過程包括三個步驟：

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

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

(3)通過適當的成像條件搜索震源位置.

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

2.1 衰減補償效果分析

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

圖4 (a) 震源與環形檢波器排列； (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 逆時定位效果比較

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

表1 三層模型參數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)彈性介質, (b)黏彈性介質無補償和(c)黏彈性介質衰減補償定位結果; 井中接收情況下(d)彈性介質, (e)黏彈性介質無補償和(f)黏彈性介質衰減補償定位結果; (g)地面接收定位結果單道數據對比，虛線為三種方法定位震源位置; (h)井中接收定位結果單道數據對比，虛線為三種方法定位震源位置; (i)地面與井中記錄衰減補償聯合定位，虛線交點為真實震源位置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 成像算子優化

圖7a為Marmousi模型，模型大小為3000 m×3000 m，網格間距10 m，將震源設置在(1500 m，1500 m)處，在地表采集地震信號，地表檢波器道間距100 m，在震源定位過程中需要考慮各種干擾，如采集信號中存在的噪聲、地下介質參數的不確定性等等.將Marmousi模型中的縱波速度、密度參數進行平滑(如圖7b、c)，設介質泊松比為0.25，計算得到橫波速度，縱波品質因子設為縱波速度的1/40，橫波品質因子設為縱波品質因子的0.83倍.在地震記錄中加入隨機噪聲，進行自相關、互相關及結合兩種方法的優化互相關算法進行震源定位.地震記錄加入強噪聲后如圖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 含噪的地表接收數據(a) 水平速度分量； (b) 垂直速度分量.Fig.8 Noisy surface data(a) Horizontal velocity component; (b) Vertical velocity component.

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

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

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

圖10 (a) 三組檢波器穿插陣列成像； (b) 兩組檢波器常規互相關成像； (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計算量情況下，常規分組成像定位結果向左方偏離，而穿插陣列分組能夠得到與原三個分組接近的定位精度.穿插陣列互相關算子能夠在噪聲壓制效果和降低計算成本上獲得更好的效果.改變模型參數平滑尺度(圖11a、b)，進行定位的結果如圖11c、d，模型參數的平滑會導致能量聚焦性減弱.在速度整體變化時，由于本方法基于波場能量進行定位，無法克服速度變化帶來的干擾，可以結合地表地震數據通過全波形反演以獲得更好的速度模型.

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

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

3.1 井中接收數據裂縫成像

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

表2 裂縫模型參數Table 2 Fracture model parameters

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

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

圖15 (a) 井中vx分量記錄； (b) 透射波記錄； (c) 散射波及層間反射波記錄； (d) 裂縫成像結果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 地面接收數據裂縫成像

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

表3 裂縫模型參數Table 3 Fracture model parameters

圖16 (a) 裂縫模型，其中“*”為震源，為檢波器； (b) 平滑模型； (c) 震源定位結果，其中虛線交點為真實震源位置，“*”為定位結果位置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中“*”為定位結果，虛線交點為真實震源位置，可以看到定位結果準確.

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

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

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

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

圖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 結論

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

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

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

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

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

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

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

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

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