摘要真实地下介质具有黏弹性，地震波在传播过程中会发生耗散与频散.忽视黏弹性介质的吸收衰减效应，逆时延拓过程中地震波将会出现振幅减弱、相位失真等现象，无法准确定位震源真实位置，因此需要对黏弹性介质中传播的波场进行衰减补偿，并通过采用合适的成像算子对微地震震源进行定位与裂缝成像.本文基于耗散与频散解耦的分数阶黏弹性波动方程模拟波场，采用lowrank分解近似混合域算子，分离衰减相关项并反转耗散项符号，并在补偿的衰减项波场的波数域中进行低通滤波，压制噪声的影响;使用优化后的成像算子进行微地震震源定位，并通过分离散射波场，对散射波进行逆时反传寻找裂缝.数值实验证明，本文方法通过lowrank近似有效提高了计算效率，衰减补偿算子在滤波器约束下能够稳定地补偿反向延拓的波场，优化后的成像算子能够在压制随机噪声的同时进一步提高计算效率和定位分辨率.
关键词黏弹性;分数阶;Lowrank;衰减补偿;微地震;逆时定位
0引言
微地震监测中震源定位是其核心问题(MaxwellandUrbancic，2001)，在弹性介质中通过波场逆时反传微地震信号能够有效确定震源位置(Steineretal.，2008;Larmatetal.，2009)，而对于黏弹性介质如果采用常规逆时反传，由于地下介质具有黏弹性会发生吸收衰减作用，引起地震波振幅减弱、相位失真，无法确定震源真实位置，因此需要对波场进行衰减补偿，采用逆时定位方法获取真实的震源位置.
黏弹性介质中微地震震源逆时定位包括两个步骤：在时间上反传检波器接收到的地震数据与补偿黏弹性介质对地震波的吸收衰减作用.黏弹性介质的速度频散和能量耗散可以采用品质因子Q描述，Kjartansson(1979)提出了常Q模型，以参考速度、衰减因子和参考频率为参数，衰减因子是频率的线性函数.进行波场逆时反传的理论基础是地震波正演模拟，近年来众多学者对黏声/黏弹性介质正演模拟的算法进行了研究.频率域数值模拟基于单频率对空间网格点求解，各频率独立计算不存在误差累积，可以方便的进行并行计算，但是在大型3D模型中该算法无法适应消息传递接口(MessagePassingInterface，MPI)并行算法及图形处理器(GraphicsProcessingUnit，GPU)加速导致效率降低;时间域数值模拟耗费内存小，操作灵活，但是存在着时间循环累积的误差.Carcione等(Carcioneetal.，2002;Carcione，2009)通过GL(Grünwald-Letnikov)近似求解时间分数阶导数，实现了黏声及黏弹性介质波动方程数值模拟，但该方法在求解时间分数阶导数时，需要保存每一时刻的应力、应变值，计算效率低，且内存需求较高.杨仁虎等(2009)提出采用拉梅差异矩阵简化方程，通过该方法模拟黏弹性介质波场具有比GL近似更高的效率.Carcione(2010)提出基于时间为二阶导数、空间为分数阶的黏声波动方程，该方程将分数阶时间导数转化到分数阶空间导数上，有效减少了内存占用，解决GL近似的效率问题.黏声方程中耗散项与频散项的解耦促进了衰减补偿逆时延拓的发展，Treeby和Cox(2010)基于分数阶黏声波动方程，将耗散项与频散项解耦，提高了计算效率，解耦的耗散项与频散项也可以用于波场反传中地震波的衰减补偿计算，有利于震源定位和偏移成像.Zhu和Carcione(2014)利用空间分数阶导数将黏弹性介质波动方程耗散与频散解耦.吴玉等(2017)将解耦的分数阶拉普拉斯算子用于黏声介质逆时偏移，改善了深层构造的成像精度.对于波数-空间域的混合域算子求解时，计算量非常大，Fomel等(2013)提出lowrank近似混合域算子，减少傅里叶变换次数提高计算效率.Song等(2013)结合lowrank频散低和有限差分法计算效率高的优点，提出了lowrank有限差分法解决声波正演问题.Sun等(2015)利用lowrank近似模拟黏声介质波场，实现了黏声介质衰减补偿逆时偏移成像;袁雨欣等(2018)采用lowrank算法求解解耦的弹性波方程，有效地减少了波场模拟数值频散，提高模拟精度.黄金强和李振春(2019)在各向异性介质中应用lowrank近似进行正演及逆时偏移，提高计算效率.近年来，为解决分数阶拉普拉斯算子计算复杂介质时成本高的问题，众多学者不断改进算法进行黏弹介质的数值模拟(Chenetal.，2016;Yaoetal.，2017;Wangetal.，2018;Guoetal.，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的黏弹性波动方程，通过lowrank近似混合域算子提高计算效率，推导衰减补偿的黏弹性介质波动方程及优化成像算子，提高震源逆时定位的精度和计算效率.在含裂缝介质中，利用模拟记录与真实记录的差异性分离散射波，基于衰减补偿和分组互相关成像算子对裂缝进行成像.
1.2Lowrank近似
在通过式(7)求解拉普拉斯算子时，从理论上讲，对于每一个c0和Q，需要为方程计算执行一次反傅里叶变换(InverseFastFourierTransform，IFFT).该方法称为逐点傅里叶变换(FastFourierTransform，FFT)方案.逐点方案是处理混合域算子的有效数值方法，但是其计算成本太大.
2微地震衰减补偿逆时定位
2.1衰减补偿效果分析
设置二维介质中震源、接收点分布如图4a，对2号环状检波器组接收到的地震记录进行逆时延拓，得到近震源1号检波器处波场如图4b，在0.12s之前的噪声是反传中互相串扰的信号以及前文提到的随机噪声，可通过成像条件进行压制.放大0.12～0.24s部分如图4c，可以看到真实波场(实线VE)与反传波场(VE-VE)存在一定能量损失，其原因是接收点不能捕获所有能量，而且为了压制噪声需要对信号进行一定程度低通滤波，也会造成少量能量损失.图中未补偿波场(VE-EL)明显出现延滞，而进行衰减补偿后波场传播时间已经校正.通过衰减补偿的lowrank波动方程进行波场延拓可以在正确位置成像，修正黏弹性介质对地震波的吸收衰减作用.
2.2逆时定位效果比较
设置一个简单的三层模型，大小为2000m×2000m，网格间距10m，参数见表1，在(1000m，1300m)处正应力分量加载主频30Hz雷克子波，分别在地表和700m处一1500m深井记录地震数据.得到多分量地震记录后(如图5所示)，可以进行单分量反传定位或多分量反传定位，地表地震记录与井中地震记录联合反传定位，也可进行分频段多尺度反传(本文默认进行全频段记录反传)等各种延拓方法.首先在黏弹性介质中分别进行弹性反传(VE-EL，图6a、d)、无补偿衰减反传(VE-QE，图6b、e)和衰减补偿反传(VE-VE，图6c、f)，并通过自相关成像算子进行震源定位，地面定位结果提取地表1000m处应力垂直分量如图6g，通过衰减补偿定位震源深度为1300～1310m，无补偿时定位在1260～1280m，使用弹性波动方程定位在1250～1280m处.无补偿波能量分布在地表到震源之间，分辨率较低，弹性波能量聚集性弱于衰减补偿波场，通过衰减补偿后波场不但定位误差最小，且能量聚集更集中便于识别;井中地震数据提取井中1300m处数据应力水平分量图6h，震源定位效果与地面定位效果基本相同，但是在等深度处400m外出现另一峰值，这是由于探测井的位置分布及井深的局限性，导致有效信息空间采样不足，残余能量聚集到井对称位置形成假震源，由于其并非真实震源，故波场能量弱于真实震源，在复杂介质中假震源聚焦能力会进一步减弱，且补偿后假震源能量更小有助于识别.通过染色算法(ChenandJia，2014)或地面记录与井中记录联合反传(图6i)等方法可以规避假震源的影响.
2.3成像算子优化
图7a为Marmousi模型，模型大小为3000m×3000m，网格间距10m，将震源设置在(1500m，1500m)处，在地表采集地震信号，地表检波器道间距100m，在震源定位过程中需要考虑各种干扰，如采集信号中存在的噪声、地下介质参数的不确定性等等.将Marmousi模型中的纵波速度、密度参数进行平滑(如图7b、c)，设介质泊松比为0.25，计算得到横波速度，纵波品质因子设为纵波速度的1/40，横波品质因子设为纵波品质因子的0.83倍.在地震记录中加入随机噪声，进行自相关、互相关及结合两种方法的优化互相关算法进行震源定位.地震记录加入强噪声后如图8，含噪信号信噪比约为-18dB，有效信号几乎被噪声淹没.
成像算子中，最大值成像法(McMechan，1982)是在逆时反传过程中不断检索波场能量最大值，以最终的最大值位置为震源空间位置，由于多次检索导致计算效率较低;自相关成像算子(Artmanetal.，2010)在波场延拓过程中进行零延时自相关，对背景噪声有一定的压制能力，但是对记录中的随机噪声无法进行有效压制(图9a);互相关成像算子(NakataandBeroza，2016)利用噪声之间相干性较低的特点，将各检波点记录单独进行逆时延拓并进行互相关，压制噪声能力较强，而随检波器的增加计算成本呈线性增长，计算效率和运行内存要额外消耗数十倍至上百倍之多，为节约计算成本一般按位置分组逆时延拓.图9b为将检波器分为等段距离的三组后进行互相关得到的定位图像，虽然压制噪声能力较自相关算法更强，但是其能量最大值位于模型上方一断层界面处.图9c为在等距三组检波器分组的基础上针对成像算子进行优化，在互相关算法的基础上进行一次自相关，结合自相关成像算子对背景噪声的压制和互相关成像算子对随机噪声压制的优点，加强成像算子的聚焦能力(Tangetal.，2020).
相关期刊推荐;《地球物理学报》创刊于1948年，是中国地球物理学会、中国科学院地质与地球物理研究所联合主办的有关地球物理科学的综合性学术刊物。主要刊载固体地球物理、应用地球物理、地磁和空间物理、大气和海洋地球物理，以及与地球物理密切相关的交叉学科研究成果的高质量论文。
对于互相关算子中检波器的分组，传统检波器阵列为[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.
对比图10a和图9c，同样计算量下，穿插阵列分组能够更好地压制噪声，而在减少1/3计算量情况下，常规分组成像定位结果向左方偏离，而穿插阵列分组能够得到与原三个分组接近的定位精度.穿插阵列互相关算子能够在噪声压制效果和降低计算成本上获得更好的效果.改变模型参数平滑尺度(图11a、b)，进行定位的结果如图11c、d，模型参数的平滑会导致能量聚焦性减弱.在速度整体变化时，由于本方法基于波场能量进行定位，无法克服速度变化带来的干扰，可以结合地表地震数据通过全波形反演以获得更好的速度模型.
2.4三维模型测试
通过三维overthrust模型进行定位测试，模型大小1500m×1500m×1200m，模型纵波速度参数如图12a所示，在地表设置四条测线进行接收(图12b)，爆炸源设置在测线交点下方深960m处.图12d为地表接收的592道垂直分量记录，在地表记录中加入信噪比为1.25dB噪声，并通过平滑后的模型(图12c)进行定位.定位结果如图13，其中图13a为通过衰减补偿定位的结果，图13b是未补偿定位结果，图13c为在真实震源处抽取的垂直单道归一化数据对比，虚线位置为真实震源深度.从图13a、b可以看到衰减补偿能量聚焦在震源附近，而未补偿定位结果聚焦效果较差，无法辨识震源;图13c中补偿后震源定位结果准确，而未进行补偿的定位结果其能量集中在地表低速带附近，震源处能量峰值位置也偏浅，存在一定误差.如果提高平滑度，即模型参数不准确时，通过衰减补偿也可能无法定位震源真实位置，因此获得准确的初始速度模型对震源定位的效果影响较高.本文方法也可以用于各向异性黏弹性介质，对于各向异性黏弹性模型的衰减补偿逆时定位结果见附录A.
3裂缝成像
在水力压裂过程中，地下存在裂缝时，根据惠更斯原理，地震波经过裂缝时会产生透射波和散射波，散射波可以视为以散射点(即裂缝)为震源发出的子波，如果能将散射波与透射波分离，对散射波和透射波进行互相关反传，它们将会在裂缝处聚焦，实现被动源裂缝成像.
3.1井中接收数据裂缝成像
设置一个层状模型对裂缝成像进行模拟，纵波速度模型如图14a，模型分为三层，大小为600m×600m，网格间距2m，在第二层存在两条裂缝，模型参数如表2，在(100m，2800m)处设置微地震震源，并于地表接收地震记录，其中水平速度分量如图15a.对图14a中的模型参数进行平滑得到的纵波速度模型如图14b所示.
通过检索每一道地震记录的走时，将初至波拉平处理;然后对地震记录进行频谱分析，通过主频与时间采样间隔估算地震波波长;再按照波长切割分离初至透射波和散射波，并通过前文介绍衰减补偿逆时反传方法进行裂缝成像，得到成像结果如图15d，井中接收数据逆时波场能量聚集在裂缝的边缘位置，误差较小.
3.2地面接收数据裂缝成像
在地表进行裂缝定位时，受层间多次波的影响，不能简单地切除初至波进行散射波分离，本文通过模拟新的波场并进行处理达到分离散射波的目的.设置纵波速度模型如图16a，模型分为四层，大小为500m×250m，网格间距1m，在第三层存在四条裂缝，模型参数如表3，在(250m，240m)处设置微地震震源，并于地表接收地震记录，其中垂直速度分量如图17a.对图16a中的速度模型进行平滑得到的速度模型如图16b所示.
(1)首先通过衰减补偿逆时定位，寻找微地震震源的空间位置，图16c中“*”为定位结果，虚线交点为真实震源位置，可以看到定位结果准确.
(2)在定位震源位置处激发震源，对图16b模型进行二次波场模拟，得到模拟记录垂直速度分离如图17b.由于已知参数不够精准，模拟直达波与真实直达波会存在走时和振幅的偏差.
(3)记录检波器每一道的最大值，将所有记录按道进行归一化处理.通过检索真实记录直达波走时，将模拟记录的直达波校正至真实直达波位置处，然后将两个记录进行相减，并将新得到的记录进行逆归一化，每道进行同样处理后得到散射波如图17c.
(4)对分离后的散射波记录与处理后的模拟波记录进行互相关反传，并切除检测到的震源能量团后成像如图17e，其中右三处裂缝位置成像能量较强，左侧裂缝由于距震源最远，成像能量较弱.
另外应用Zhu(2019)提出的中值滤波方法进行测试，得到散射波如图17d，成像结果如图17f，其结果与本文方法结果相似，左侧裂缝能量较弱.
4结论
针对黏弹性介质震源定位及裂缝成像中计算效率及准确度的问题，本文基于Kjartansson常Q模型及Carcione的黏弹性介质波动方程，引入lowrank算法进行计算加速，并对成像算子进行了改进，实现衰减补偿震源定位及裂缝成像，得到以下结论和认识：
(1)通过分数阶拉普拉斯算子将能量耗散和相位频散明确分离，在逆时反传时通过反转耗散算子的符号并使保持频散项符号不变来完成衰减补偿.补偿后的波场能够恢复能量并校正能量聚集位置，定位震源真实位置.
(2)针对计算混合域算子时逐点傅里叶变换计算成本较高的问题，通过lowrank分解将混合域算子分解为三组矩阵，采用巴特沃兹低通滤波器稳定衰减补偿波场，避免了高频噪声被补偿引起的高能量波场掩盖有效信号，在保证波场模拟精度的情况下降低了计算成本.
(3)针对逆时定位过程中存在的随机噪声影响，优化成像算子，在减少计算量的同时提高定位分辨率.优化后的成像算子能够在加入强随机噪声的地震记录中高效率高精度定位震源位置.
(4)对于水力压裂过程中介质中存在的裂缝，通过分离地震记录中的透射波和散射波并进行互相关反传，能够对裂缝进行成像，但实际上常规方法对散射波的分离效果仍不理想，新近发展的深度学习方法可能是一个解决方向.
由于采用常规的波场能量作为成像值，在速度存在整体偏离的情况下震源位置存在误差.对于速度整体的偏离引起的误差可以与全波形反演方法结合，本文采用的成像算子也可以移植到能量范数、反褶积算法等其他算法上.震源定位问题涉及初始模型参数、各向异性、黏弹性以及震源机制等多方面，本文主要解决黏弹性介质带来的地震波衰减的影响，并未对其他因素详细讨论，涉及多参数影响下的震源定位问题，或许可以借助全波形反演的思想进行解决，尚需进一步研究.
附录A三维黏弹性各向异性介质震源定位
本文讨论了二维各向同性黏弹性介质的衰减补偿逆时定位，该方法在各向异性介质中同样适用.这里设计了一个黏弹性VTI介质进行测试，对离散采样的overthrust模型(如图A1a)添加速度各向异性，设置Thomsen参数ε=0.2，δ=0.1，γ=0.15，震源为加载在中心点深900m处的爆炸源，地表设置四条测线，利用平滑后的模型进行逆时定位，成像算子采用本文设计的优化互相关算子.图A1b为检波器及震源位置，图A1c为地表592道垂直分量记录，图A1d为震源定位结果剖面，图A1e为抽取震源位置垂直单道归一化曲线，可以看到，在模型中添加速度各向异性并不影响本文方法的定位效果，在模型参数较准确时，本文方法能够准确定位震源.——论文作者：唐杰，刘英昌，李聪，孙成禹
转载请注明来自：http://www.lunwencheng.com/lunwen/lig/21005.html
