基于成像算子优化的微震逆时定位方法

唐 杰， 李 聪， 孙成禹， 陈学国

(1.中国石油大学(华东)地球科学与技术学院，山东青岛 266580；2.中国石油化工股份有限公司胜利油田分公司勘探开发研究院，山东东营 257015)

水力压裂能够改善储层的孔渗特征，对于致密油气的增产是十分有效的[1-2]。在压裂过程中岩石破裂会诱发微地震，通过微地震的监测可以有效评估压裂的效果[3-4]。在微地震监测中，震源定位作为关键问题，目前主要有基于走时和波形的定位方法[5]。Geiger最早通过到时信息确定了天然地震的位置[6]，后续发展的双差定位方法、各种走时计算方法都是在此基础进行的优化[7-10]。这类基于走时信息微震定位方法需要高信噪比的微震记录来精确提取初至波的走时信息。微震逆时定位方法通过波场的反向延拓，借助不同的成像条件，使波场在震源位置聚焦，由于在计算过程中使用的是全部波场的信息，不用对走时信息进行拾取，适合处理低信噪比的地面观测资料[11-16]。Sava 等[17-18]提出了一种干涉成像条件可以有效压制波场串扰，提高成像精度。Li等[19]开发了一种基于定理表示的微震定位方法，使在反向传播过程中的波场定向聚焦到震源的位置。Douma和Snieder[20]使用反褶积成像条件，提高了微震逆时成像的空间分辨率。Nakata和Beroza[21]提出了互相关成像条件，得到具有较高空间分辨率的震源成像点。李萌等[22]提出以峰值信噪比和定位椭球半轴长作为定量评价各种因素对定位结果的影响。徐克彬等[23]采用L-M反演算法逐次迭代寻找能量叠加最大值点。李青峰等[24]提出分组互相关成像条件可以在满足计算效率要求的情况下,获得空间分辨率很高的震源成像结果。笔者以传统的成像条件为基础，对传统成像算子的定位结果缺陷进行分析，提出一种优化后成像条件，分别对简单模型和复杂模型进行测试。

1 微震逆时成像定位

微震逆时定位是将微震记录的波场逆时反传成像。地面各检波点接收的地震记录是由同位置的微地震震源激发得到，通过借鉴逆时偏移成像的思想，将各个检波点得到的波场记录进行反向传播，采用合适的成像条件就可以将震源位置在空间上聚焦出来。微震逆时定位的两大关键步骤为波场反向延拓和成像条件的选择。

1.1 波场逆时延拓方法

逆时延拓可以看作是波场正演模拟的逆过程。假设震源位置和检波点位置分别为Xs和Xr，波场正演过程可以表示为

U(Xr,t)=F-1{S(Xs,ω)G(Xr,Xs,ω)}.

(1)

式中,G和S分别为格林函数和震源函数;t和ω分别为传播时间和频率;U为检波器记录波场;F-1为傅里叶逆变换。

波场的逆时延拓便可表示为

(2)

式中,UR为逆时延拓波场;X为地下空间位置,2维情况下为(x,z);U为检波点ri接收到的波场,作为波场逆时延拓的震源,i表示第i个检波点。

研究采用的是时间二阶、空间八阶的有限差分方法进行波场的正演和反传过程。

在震源处的二维声波方程可以表示为

(3)

式中,Us(X,T)为从震源处激发产生的波场,即正演波场;S(Xs,T)为震源项;T为正演过程传播时间;X为空间位置(x,z)。

在波场反传过程中,将震源项替换为检波器处的波场,时间进行反向传播,可获得声波逆时延拓波场[13],公式为

(4)

式中,UR为逆时延拓产生的波场;U为检波器处的波场;Xr为检波器的位置;v为介质空间的速度;T为微震记录总时间长度;t为逆推波场传播时间。

1.2 成像条件

1.2.1 最大振幅成像算子

最大振幅成像算子是最常用的逆时定位成像条件,可以表示为

I(x,z)=maxt=Tmax|UR(x,z,t)|.

(5)

式中,UR(x,z,t)为地下空间位置(x,z)和传播时间t对应的逆时延拓波场;I(x,z)为成像剖面中的成像值。

提取出逆时延拓过程中所有时刻的最大振幅作为震源位置的成像值。

1.2.2 自相关成像算子

自相关成像将各检波点微地震记录作为震源函数,进行波场逆时延拓计算。这种延拓方法对成像结果中背景的噪声有一定的压制效果,在波场逆时延拓过程中的自相关成像算子可表示为

(6)

式中,Tmax为记录微地震波场的最大时间;I(x,z)为成像剖面中成像值。

1.2.3 互相关成像算子

地震记录中的噪声大多为随机分布,噪声与噪声、噪声与信号之间的相干性很低,而信号和信号之间的相干性较高。互相关成像条件把各个检波点作为相互独立的震源,分别进行波场逆时延拓,波场在回到发震位置处的相干程度最大,因此在震源处可以得到最大的波场互相关值,同时又能压制波场中互不相干的噪声,成像结果的分辨率较高。在波场逆时延拓过程中应用的互相关成像条件[20]可表示为

(7)

式中,URi(x,z,t)为第i个检波点在空间位置(x,z)和传播时间t的逆时延拓波场;N为所用的检波点的总个数。

小学语文的课文通常都包含着作者浓烈的思想感情，这份感情是显而易见的，学生们若是朗读课文，便更能体会到作者的感情，甚至可以达到感同身受身临其境的境界。朗读也可以让学生熟悉作者在写作中所运用的写作技巧，从而提高学生的语文素养，学生可以运用到平日写作中。可见，朗读对小学生学习小学语文有着直接且效果明显的帮助。

1.2.4 优化后的成像算子

最大振幅成像算子在成像过程中直接对反传波场的振幅进行叠加,低信噪比的微震记录中含有较多噪音,直接进行波场的叠加会导致大量噪音也会叠加到成像剖面中,导致成像的分辨率降低,对定位结果有一定程度影响。自相关成像算子在波场逆时延拓过程中也会对波场延拓过程中震源外的干扰波场进行自相关,会对定位成像结果产生一定程度的干扰。互相关成像算子虽然可以对互不相干的随机噪音进行压制,但是当噪声水平较高时,需要更多的检波点来压制噪声,并且压制自身波场延拓产生的误差能力较弱。

针对传统的成像算子存在的问题,提出了优化的成像算子,在基于互相关成像条件的基础上对其进行一次自相关,逆时定位过程中的成像条件可表示为

(8)

式中，Tmax为记录微地震波场的最大时间。

自相关成像算子和互相关成像算子在波场逆时定位中对提高定位成像效果的作用是不同的，不同检波器的逆时延拓波场最终在震源处会产生一个相干加强的波场最大值，通过自相关成像算子可以凸显出波场中的有效信息，一定程度上压制背景噪音，互相关成像算子主要根据随机噪音互不相干的特点，通过互相关算法对成像后道间互不相干的噪声进行压制，同时也可以使震源处的相干波场最大化。优化后的成像算子充分结合了这两种成像算子的优点，最大化地压制了波场中的背景以及随机噪音，同时由于自相关算子对噪声水平有一定的压制效果，可以互补互相关成像算子需要高密度检波点的缺点，在保证检波器数量的前提下可以使定位剖面中震源处的有效信息得到高质量的显示。

1.3 逆时定位结果评价标准

由于定位图像中的能量会在震源处聚焦，通过提出采用统计学中峰度的概念来对定位聚焦效果进行定量评价。峰度[25]是指通过与正态分布进行比较来描述总体数据取值分布形态陡缓的程度，峰度越大说明定位结果更加集中，定位成像效果更加清晰，与之相反，峰度越小表明定位结果分布更为分散，成像质量不佳，分辨率较低。峰度计算公式为

(9)

2 模型测试

2.1 简单层状模型测试

模型选用三层层状模型，大小为2 000 m×2 000 m，网格大小为5 m×5 m，检波器设在地表，间隔为5 m，共400个检波器，模拟过程中采用声波有限差分法进行正演和逆时波场反传。模型如图1(a)所示，x和z分别表示水平方向和深度方向，“+”表示单源位置，微震正演模拟结果如图1(b)所示。

图1 层状介质模型及微震数据Fig.1 Layered media model and micro-seismic data

图2 层状模型4种成像算子定位效果Fig.2 Positioning effect of four imaging operators in layered model

为了对4种成像算子的成像效果有定量的评价，将真实震源处水平和深度方向的峰度给出，图3为真实震源处各种成像条件在水平和深度方向上的峰度。从图3中可以看出，优化后的成像算子在深度和水平方向上的峰度都要高于其他3种成像算子，说明采用优化后的成像算子在定位成像剖面中的震源点聚焦效果更好，成像质量相对较好，成像精度也相对较高。

图3 简单模型在震源处水平和深度方向的峰度Fig.3 Kurtosis in horizontal and depth directions of simple model at source

针对简单层状模型的测试，优化后的成像算子可以较好地给出震源的位置，成像分辨率和精度也得到了提高。

2.2 复杂模型测试

2.2.1 复杂模型测试

野外的地下地质情况复杂多变，为了模拟野外实际情况，对优化后的成像算子进行复杂模型的测试。本文中选取的复杂模型为部分Marmousi模型，模型大小为2 000 m×2 000 m，网格大小为5 m×5 m，检波器设在地表，间隔为5 m，共400个检波器。模型如图4所示，图中“+”表示单源位置。

图4 复杂介质模型及微震数据Fig.4 Complex medium model and micro-seismic data

图5为复杂模型采用4种成像算子得到的定位效果。对比分析成像结果，4种成像算子在定位剖面的上半部分由于模型的复杂性都出现了一定程度的干扰，但采用优化后的成像算子能量依然较好地在真实源位置得以集中，相比较其他3种成像算子的成像效果，定位剖面较为清晰，效果图上部干扰在一定程度上得以压制。

图5 复杂模型4种成像算子定位效果Fig.5 Positioning effect of four imaging operators in complex model directions of complex accurate model at source

图6为复杂模型采用4种成像算子定位后在真实震源处的水平方向上和深度方向上的峰度值。对比分析4种成像算子的峰度值，可以看出，采用优化后的成像算子的峰度值在水平和深度方向上的峰度值都高于其他3种传统成像算子，说明成像定位能量聚焦较好，在深度和水平方向上的成像分辨率也相较于传统成像算子在一定程度上得到提高，针对复杂模型有较好的成像结果。

图6 复杂精确模型在震源处水平和深度方向的峰度Fig.6 Kurtosis in horizontal and depth directions of complex accurate model at source

2.2.2 复杂模型含噪数据测试

对复杂模型的微震数据加入信噪比为1的随机噪音并对成像算子进行抗噪性测试，图7为加入随机噪音后的正演数据。图8为复杂模型含噪数据采用4种成像算子的定位效果。

图8 复杂模型含噪数据4种成像算子定位效果Fig.8 Positioning effect of four imaging operators in complex model with noisy data

对比分析4种成像结果，可以看出复杂模型的含噪数据对定位影响的效果比较大，特别针对传统成像算子，定位结果几乎淹没在成像噪音的干扰中，优化后的成像算子成像剖面虽然受到一定程度的影响，但是震源位置依然可以清晰分辨出，震源处的能量得到较好的聚焦，定位成像效果图震源位置分辨率也较高。

图9为在真实震源处提取水平方向和深度方向上的4种成像算子的峰度值，可以看出复杂模型含噪数据对深度方向上的自相关和互相关成像算子的定位效果产生一定程度的影响，但是采用优化后的成像算子的定位结果依然可以保持较高的峰度值，说明定位剖面的能量聚焦较好，成像剖面的分辨率较高，定位结果的精度相较于传统成像算子有一定程度的提高，证明优化后的成像算子具有较好的抗噪性。

图9 复杂模型含噪数据在震源处水平和深度方向的峰度Fig.9 Kurtosis in horizontal and depth directions of noisy data in complex model at source

2.2.3 平滑复杂模型测试

在运用实际数据进行逆时定位过程中所采用的地下介质速度通常都是不精确的，因此为了模拟真实情况下的速度不确定性，对复杂的速度模型进行平滑处理，波场反传过程中采用平滑后的速度模型，复杂平滑速度模型如图10所示。

图11为针对平滑后的复杂模型采用4种成像算子得到的定位效果。从图11中可以看出平滑后的复杂模型对几种成像算子的定位效果影响都不大，但是相比其他3种成像算子，采用优化后的成像算子定位精度更高，成像剖面中干扰压制更好，能量更加集中。

图10 平滑复杂介质模型Fig.10 Smooth complex media model

图12为针对平滑复杂模型采用4种成像算子在真实震源处提取的水平方向上和深度方向上的峰度。可以看出优化后的成像算子在震源处水平和深度方向上相比于传统成像算子都能保持较好的聚焦，震源处横向和纵向的分辨率都相对较高，证明优化后的成像算子对速度模型的不确定性有较好的适应性。

图11 平滑复杂模型4种成像算子定位效果Fig.11 Positioning effect of four imaging operators in smooth complex media model

图12 平滑复杂模型在震源处水平和深度方向的峰度Fig.12 Kurtosis in horizontal and depth directions of smooth complex media model at source

通过对复杂模型进行的一系列测试表明，采用优化后的成像算子能较好地适应不精确的模型以及低信噪比的微震数据，对比传统成像算子，本文中提出的优化后的成像算子既能获得准确位置，在精度上又可以得到提高，并且也有广泛的适用性。

2.3 检波器分布及数量的测试

2.3.1 检波器分布测试

野外观测系统的分布与理论设计的观测系统通常有较大的差异，因为在野外观测时通常会受到地形、地貌等一系列因素的影响，导致观测系统中的检波器排列不均，针对检波器的分布问题对测试优化后的逆时定位成像算子进行测试。选用3种排列分布的检波器进行测试分析，检波器分布如图13所示。

图14为采用优化后的成像算子对3种分布的检波器的逆时定位效果。对比全分布的检波器定位效果，部分分布的检波器的定位效果受到了一定程度上的影响，并且采用不同分布方式的定位效果受到的影响程度还有差异。图14(b)和(c)受到的影响程度较大，定位剖面上部受到了很强的干扰，这可能是模型的复杂性所导致的，但是无论采用哪种分布，在震源处都有很好的能量聚焦，震源的位置都可以确定。

图13 3种检波器分布范围示意图Fig.13 Schematic diagram of three kinds of receivers distribution range

图14 3种分布范围的检波器定位效果Fig.14 Positioning effect map of three kinds of receivers distribution range

2.3.2 检波器数量测试

检波器的数量通常决定着成像的精度，互相关成像算子虽然能较好地压制不相干的随机噪音，但是需要更多的检波点进行压制，优化后的成像算子在基于互相关成像算子上进行改进，其中一个优点就是可以减少检波器的使用数量。

检波器数量的测试选用图4(a)部分Marmousi复杂模型进行测试。为了测试结果有定量的评价，通过选取成像最大值位置与真实位置的误差距离来衡量评价定位结果,

E=|(x,z)Imax-(x,z)t|.

(10)

式中，E为绝对误差距离；(x,z)Imax为最大成像值的位置；(x,z)t为震源的真实位置。

图15是针对4种成像算子分别采用不同数目的检波器进行逆时定位得到的位置与真实位置的绝对误差距离。通过对比4种成像算子与检波器数量的关系可以得到，当检波器数量较多时定位得到的震源位置都比较准确，误差较小甚至没有误差；随着检波器数目的减少，采用最大振幅成像算子最先开始产生较大误差，采用优化后的成像算子在检波器减少到一定程度上依然可以保持较好的成像精度，同时在保持成像精度的前提下所需要的检波器数目最少。所以采用优化后的成像算子可以较好地适应检波器数量较少的情况，并且相比较传统成像算子有更好的成像结果。

图15 4种成像算子不同检波器数量的绝对误差距离Fig.15 Absolute error distance of four imaging operators with different numbers of detectors

3 结 论

(1)优化后的成像算子相比较于传统成像算子在定位能量聚焦程度和成像分辨率上都得到了提升，对无源处的干扰特别是浅部干扰有较好的压制效果，对含速度误差和含噪微震数据都有较好的成像结果。

(2)采用优化后的成像算子可以较好地适应检波器不同分布方式和数量较少的情况，能够获取准确的震源位置。