基于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