基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法与流程

本发明属于地震勘探技术领域，涉及一种反假频射线束形成偏移成像方法，尤其是一种基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法。

背景技术

地震信号处理所面对的是离散的多维数据，尤其是三维情况下，假频是必须面对的问题之一。三维叠前地震数据本质上为五维信号，包括时间，二维炮点位置坐标和二维检波点坐标等五个维度。同样，三维偏移成像结果也可为五维数据，即空间位置坐标和二维水平坐标以及地下层位的反射角和方位角。因此，三维叠前偏移成像是多维数据到多维数据的映射，而且每个维度上均为离散量，会涉及多个采样间隔。首先是观测系统间隔，这其中又涉及到炮点网格和检波器网格，分别在二维空间上采样；其次是地震成像网格，即输出地震图像三维空间位置采样以及反射角和方位角采样间隔。理论上，每个维度上的离散化均需分别满足采样定理。然而，由于受采集成本或者施工限制影响以及综合考虑数据存储量、计算量等因素，无论是数据采集还是后续数据处理环节，均在尽量大的网格上进行操作，不能够总是满足采样定理的要求，即会引入假频。一般称此类由采样引起的假频为数据假频。内插是压制数据假频的有效方法(刘喜武等，2004；刘财等，2013)，可以解决观测系统不规则及由偏移引起的部分假频问题。但插值重建涉及较大计算以及存储量(刘礼农等，2006)，数据量的增大直接降低后续数据处理尤其是偏移成像的计算效率。

数据假频是引起偏移成像中假频噪音的一个原因，另一原因是由叠前偏移操作引起的，称为算子假频(zhouetal.,2013；nimsailaetal.,2015)。地面接收地震数据是地下介质反射或者散射波场的采样，各种积分类偏移即可视为在叠前五维地震数据中有选择性的“采样累加”，如kirchhoff偏移(schneier,1978)和各类射线束偏移(grayetal.,2009)。但对偏移时的地震道采样间隔，当偏移算子斜率增大到一定程度后，并不能满足采样定理从而导致假频噪声。算子反假频在kirchhoff偏移中得到了较为充分的发展。最简单高效的假频压制方法即是限制偏移孔径(王华忠等，2010)或者直接限制偏移算子的大角度分量，该策略对小角度反射层结构成像改善效果明显，但并不能完全消除偏移假频，而且对成像陡倾角结构损伤巨大(abma,etal.,1999)。原则上，反假频需要计算旅行时对空间的导数，该操作在偏移的最内层循环，会严重降低偏移的计算效率(zhangetal.,2001)。为了计算简单以及提高偏移效率，可以用常速(地表速度)来代替计算旅行时对空间的导数。虽然这种近似操作对高倾角结构成像的分辨率和振幅保真度有所下降(zhouetal.,2013；gray,2013)，但由于其实用高效，仍然在kirchhoff偏移反假频中得到广泛应用。gray(1992)提出与kirchhoff偏移算子倾角相关的低通滤波反假频算法，该方法利用不同频带的低通滤波器对地震道滤波，产生多个不同带宽的地震数据，偏移时根据算子倾角选择不同数据作为输入。为了实现输入数据在偏移剖面上的平滑过渡，gray的方法对每道数据需要产生一定数量的截止频率不同的版本，增加了数据预处理计算量和存储量。claerbout(1992)和lumleyetal.(1994)利用n点三角滤波器，提出了一种实时计算的偏移算子反假频kirchhoff偏移，以增加一定计算量的代价避免了多频带数据的产生和存储。王华忠等(2010)给出了三维kirchhoff共偏移距叠前时间偏移中的反假频准则。针对上述反假频中近似带来的误差，pica(1996)、biondi(2001)以及zhouetal.(2013)提出了数据倾角自适应反假频方法，对目标倾角结构成像分辨率和振幅保真均有提升，但牺牲了非目标倾角结构成像质量。nimsailaetal.(2015)提出在频率波数域或者频率慢度域进行炮域数据反假频滤波，用来压制高频逆时偏移成像中的假频噪音。

射线束成像在计算效率、陡倾角成像、信噪比、算法稳定性以及灵活性等方面均有一定优势。该方法将邻近道数据同时向下偏移，对宽方位、宽频带和高密度地震数据采集偏移成像优势明显。尤其是共偏移距射线束偏移成像，地震信号线性特征强，利于射线束形成以及输出偏移距道集。波束形成是射线束偏移中关键步骤，利用各种局域化分解技术将地震数据分解为不同方向的局部平面波(波束)，提取出来的局部平面波是叠前地震数据的一类特征波场，其结果直接影响偏移效率和质量。所以有必要针对共偏移距射线束偏移方法，提出一种基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法。

技术实现要素：

本发明的目的在于克服现有技术的不足，提出了一种基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法。该方法中三角反假频滤波算子直接在时域实时完成，在对地震数据进行因果和反因果积分后，可以直接与时域倾斜叠加波束形成结合完成，避免了将数据变换至频域进行滤波操作。同时，本发明中提出在反假频公式中加入权重系数，以对反假频程度进行控制。

本发明的目的是通过以下技术方案来实现的：

基于三角滤波的局部倾斜叠加波束形成偏移假频压制方法，包括以下步骤：

(1)对给定的三维叠前反射数据u(xs,ys,xr,yr,t)，对数据划分形成超道集

其中，xs,ys分别是炮点在联络线和主测线方向上的位置坐标，xr,yr分别是检波点在相应方向上的位置坐标，t是时间坐标，分别是每一道数据相对于超道集中心位置的局部相对坐标，分别是超道集炮点中心位置坐标，分别是超道集检波点中心位置坐标；

(2)对形成的超道集进行坐标变换，变换后的超道集利用局部中心点局部偏移距坐标表示为

(3)对坐标变换后的超道集进行部分正常时差校正，消除局部偏移距δhx和δhy的影响，得到部分nmo后的超道集数据

(4)对nmo后超道集进行因果和反因果积分后，进行时移叠加即可实现带反假频效果的倾斜叠加波束形成结果u(pmx,pmy,τ)。

本发明进一步的改进在于，步骤(2)中坐标变换为

δhx＝δxr-δxs，δhy＝δyr-δys

式中，δhx和δhy分别表示局部偏移距，δmx和δmy分别表示局部中心点坐标。

本发明进一步的改进在于，步骤(3)中进行部分正常时差校正采用如下公式进行

此公式表示需将超道集地震道(δhx,δhy)处t(δhx,δhy)时刻波场平移至使该地震道数据等效为由定义的偏移距接收数据；表示中心点为处的rms速度，t0表示零偏移距时间；其中坐标和是由超道集中心的源-检波器坐标转换为偏移距-中心点坐标，定义如下

本发明进一步的改进在于，步骤(4)的具体过程为：

(1)对部分正常时差校正后的超道集数据进行因果积分和反因果积分，得到

(2)对给定的中心点射线参数(pmx,pmy)，计算倾斜叠加时移量t_shift＝pmxδmx+pmyδmy；

(3)对给定的中心点射线参数(pmx,pmy)，计算反假频滤波算子时移量以及归一化因子α＝(k+1)²，其中，最大频率dt是时间采样间隔，dx和dy分别是联络线和主测线方向上的空间采样间隔，wx和wy是引入的反假频权重系数，分别控制联络线和主测线方向上反假频量；

(4)对利用偏移倾斜叠加时移量和反假频滤波算子时移量，累加至倾斜叠加结果

其中，τ是倾斜叠加的截距时间。

与现有技术相比，本发明具有以下有益的效果：

(1)本发明的三角反假频滤波算子可以直接在时间域实现，避免了将数据变换至频域进行滤波操作，从而无需额外的正反傅里叶变换。

(2)本发明的提出的反假频滤波算子的因果和反因果积分操作，可以作为预处理应用于地震数据，从而可以直接与时域倾斜叠加波束形成结合。

(3)本发明提出的反假频滤波算子引入了反假频权重系数，提高了对反假频操作的可控性，避免了三角滤波反假频操作对高频数据的过分压制。

(4)本发明提出的反假频方法是在叠前数据域完成，根据地表速度决定的慢度分量在数据域滤波，反假频后数据可进行压缩存储后重复使用，不受速度模型更新影响。

附图说明

图1是积分法偏移算子假频产生示意图。

图2是海上地震数据某超道集中源、检波点位置示意图，该超道集中共有367条数据；其中，(a)源点位置；(b)检波点位置。

图3是超道集数据及部分nmo：其中，(a)原始数据；(b)部分nmo后数据。

图4是k取3、4、5的三角滤波器频谱图。

图5是三维倾斜叠加结果，其中crossline和inline角度范围均为-65°～65°，角度间隔为5°。其中，(a)无反假频；(b)反假频权重系数wx＝1.0，wy＝1.0；(c)反假频权重系数wx＝1.5，wy＝1.5。(d)选取(a)图中部分慢度系数放大显示；(e)选取(b)图中部分慢度系数放大显示；(f)选取(c)图中部分慢度系数放大显示。

图6是某海上实际数据共偏移距射线束三维成像一条inline线结果；其中，(a)无反假频波束形成；(b)反假频波束形成(反假频系数wx＝1.0，wy＝1.0)；(c)反假频波束形成(反假频系数wx＝1.5,wy＝1.5)。

图7是海上地震数据三维成像结果(图6)选取两部分区域波形图放大显示(区域一和二位置如各图中坐标所示)；其中，(a)为所选区域一无反假频波束形成；(b)为所选区域二无反假频波束形成；(c)区域一反假频波束形成(反假频系数wx＝1.0，wy＝1.0)；(d)区域二反假频波束形成(反假频系数wx＝1.0，wy＝1.0)；(e)区域一反假频波束形成(反假频系数wx＝1.5,wy＝1.5)；(f)区域二反假频波束形成(反假频系数wx＝1.5,wy＝1.5)。

图8是共成像点偏移距道集；其中，(a)无反假频波束形成；(b)反假频波束形成(crossline和inline方向反假频权重系数wx＝1.0，wy＝1.0)；(c)反假频波束形成(crossline和inline方向反假频权重系数wx＝1.5,wy＝1.5)。

图9是三维偏移深度切片；其中，(a)深度500m；(b)深度3km；(c)深度5km。

具体实施方式

下面结合附图对本发明做进一步详细描述。

三维叠前地震数据u(xs,ys,xr,yr,t)本质上是一个五维的信号，包括炮点在crossline和inline方向上的位置坐标xs,ys，检波点在相应方向上的位置坐标xr,yr，以及时间t等五个维度。对此数据进行各种积分类偏移可视为在叠前五维地震数据中有选择性的“采样累加”。图1所示为二维叠前积分偏移示意图，从偏移算子累加曲线上可以看到，偏移算子在各个相邻道之间的等效时间采样间隔δtn,(n＝1,2,3,4,5)，不仅与地震道的采样间隔δx有关，还与偏移算子累加曲线的斜率有关。当偏移算子斜率增大到一定程度后，就会导致算子假频。针对射线束偏移，在倾斜叠加波束形成时结合三角低通滤波，可以实现偏移假频的压制。

射线束偏移将一定数量的相邻道数据同时偏移，这些同时向下传播的地震道组成的集合称为超道集(supergather)。设超道集的炮点中心位置坐标分别是超道集检波点中心位置坐标分别是则形成超道集后，三维叠前反射数据u(xs,ys,xr,yr,t)可用相对于超道集中心位置坐标的局部坐标表示

式中，分别是每一道数据相对于超道集中心位置的局部坐标。将超道集简记为对其中炮点和检波点均做倾斜叠加的时间域公式为

其中，(psx,psy)是震源的射线参数，(prx,pry)是检波器射线参数，τ表示倾斜叠加的截距时间。

对形成的超道集可以通过坐标变换，将局部坐标δxs，δys，δxr和δyr转换为局部中心点局部偏移距坐标。此坐标变换为：

δh＝δx-δx

xrs

δhy＝δyr-δys

式中，δhx和δhy分别表示局部偏移距，δmx和δmy分别表示局部中心点坐标。经过坐标变换后的超道集表示为则对应的倾斜叠加公式可写为：

其中，(phx,phy)是偏移距射线参数，(pmx,pmy)是中心点射线参数，该组参数与震源射线参数，检波器射线参数的关系如下：

phx＝prx-psx

phy＝pry-psy

在偏移距-中心点坐标系下，波束形成可针对偏移距坐标或者中心点坐标叠加，即共中心点波场或者共偏移距波场波束形成。共偏移距数据类似一个三维叠后剖面，其中反射信息连续性和相关性好，更易于倾斜叠加提取其中线性信号，因此本发明在共偏移距域实现波束形成和偏移。实行共偏移距波束形成之前，需将超道集中的局部偏移距(δhx,δhy)消除，即将所有道偏移距校正至以超道集中心定义的偏移距上，需要对超道集进行部分nmo。对超道集数据进行部分nmo校正公式为

此公式表示需将超道集地震道(δhx,δhy)处t(δhx,δhy)时刻波场平移至使该地震道数据等效为由定义的偏移距接收数据。表示中心点为处的rms速度，t0表示零偏移距时间。其中坐标和是由超道集中心的源-检波器坐标转换为偏移距-中心点坐标，定义如下

经过部分nmo校正之后的超道集已经没有了局部偏移距的差异，也即δhx＝δhy＝0。因此，坐标变换后的超道集数据可进一步简化为部分nmo校正后的超道集相应的倾斜叠加公式为

针对此共偏移距倾斜叠加公式，本发明提出利用三角低通滤波算子进行反假频，三角滤波器的z变换域表示为

其中α＝(k+1)²为归一化因子，是滤波器算子长度的一半，由最大频率决定，其中，dt是时间采样间隔，dx和dy分别是crossline和inline方向上的空间采样间隔，wx和wy分别是引入的反假频权重系数。从而偏移算子积分反假频准则为：

其中f是偏移曲线上不同道处数据累加所能满足采样定理的频率。图4为k取3,4,5时的三角滤波器频谱图，由图中可见三角滤波器可有效衰减数据中的高频成份，从而达到滤波的目的。对三角滤波器，其中每个分式可视为独立部分，三部分可以对地震道按顺序依次独立进行，而g1(z)＝1/(1-z^-1)和g2(z)＝1/(1-z)分别表示对地震道的因果和反因果积分且与k无关，因此可以作为预处理应用于地震数据。对部分nmo校正后的超道集首先进行因果积分和反因果积分，得到然后作用第三项由z变换的性质，其转换到时间域可以通过差分方程实现。从而，反假频的倾斜叠加结果为

其中，是，t_shift＝pmxδmx+pmyδmy是倾斜叠加时移量，τ是倾斜叠加的截距时间。

本发明旨在针对三维叠前地震数据的射线束偏移，提出一种基于时间域三角滤波的局部倾斜叠加波束形成偏移假频压制方法。其具体步骤如下：

(1)对给定的三维叠前反射数据u(xs,ys,xr,yr,t)，对数据划分形成超道集u(δxs,δys,δxr,δyr,t)

其中，xs,ys分别是炮点在crossline(联络线)和inline(主测线)方向上的位置坐标，xr,yr分别是检波点在相应方向上的位置坐标，t是时间坐标，分别是每一道数据相对于超道集中心位置的局部相对坐标，分别是超道集炮点中心位置坐标，分别是超道集检波点中心位置坐标。

(2)对形成的超道集进行坐标变换，变换后的超道集利用局部中心点局部偏移距坐标表示为其中坐标变换为

δhx＝δxr-δxs，δhy＝δyr-δys

式中，δhx和δhy分别表示局部偏移距，δmx和δmy分别表示局部中心点坐标。

(3)对坐标变换后的超道集进行部分正常时差校正(nmo)，消除局部偏移距δhx和δhy的影响，得到部分nmo后的超道集数据

(4)对部分nmo后的超道集数据进行因果积分和反因果积分，得到

(5)对给定的中心点射线参数(pmx,pmy)，计算倾斜叠加时移量t_shift＝pmxδmx+pmyδmy；

(6)对给定的中心点射线参数(pmx,pmy)，计算反假频滤波算子时移量以及归一化因子α＝(k+1)²，其中最大频率dt是时间采样间隔，dx和dy分别是crossline和inline方向上的空间采样间隔，wx和wy是引入的反假频权重系数，分别控制crossline和inline方向上反假频量。

(7)对利用偏移倾斜叠加时移量和反假频滤波算子时移量，累加至倾斜叠加结果

其中，τ是倾斜叠加的截距时间。

实际数据测试结果

本测试所用工区inline和crossline方向长度分别为20.75km和13.75m，数据时间采样间隔4ms，记录长度6s，最大偏移距6km。某一超道集源-检波器位置分布如图2所示，该超道集inline方向偏移距为600m，crossline偏移距0m，crossline和inline采样间隔分别为25m和12.5m。超道集大小的选择一般根据模型复杂度、射线束参数设置、目标层深度以及成像精度要求和计算时间综合考虑。选择超道集大小的一个直接有效方法即是应用不同大小超道集划分对同一块测试数据进行多次偏移。选择一条inline线作为测试线，所用数据孔径为crossline方向全数据，inline范围为该测试线两侧各4km内数据，所选择孔径范围内数据总大小为352g。综合考虑计算能力和成像质量，本次偏移超道集大小crossline和inline均为100米，角度范围均为-40°到40°，采样间隔10°，倾斜叠加后数据经压缩存储后为432g，是原始数据的1.22倍。

图3(a)和(b)分别为原始和nmo后超道集，二者在大时间处只有微小差异，但是在浅层小时间处(小于1s)图3(a)中略有倾斜的反射信号在图3(b)中被校平。

图5所示为图3中超道集数据倾斜叠加结果，其中crossline和inline角度范围均为-65°～65°，角度间隔为5°。图5(a)为无反假频倾斜叠加结果，为三维数据体。为了显示成二维图像，设置慢度坐标pmx为快维，pmy为慢维，可以看出明显的随机噪音分布在所有角度分量上。图5(b)和图5(c)为反假频权重系数全部取1.0和1.5的倾斜叠加结果，应用三角反假频后随机噪音得到明显压制。图5(c)相对5(b)在大角度能量上衰减更强，反假频后大角度倾斜叠加数据分辨率明显低于小角度。为了显示细节，图5(d)、5(e)和5(f)分别选取图5(a)、5(b)和5(c)中的部分慢度倾斜叠加结果放大显示。由于本超道集数据采样稀疏，倾斜叠加不仅有表现为随机噪音的算子假频，还有在大角度出现的数据假频(如图5(d)中椭圆标识部分)。反假频操作不仅压制了偏移假频，而且还能压制部分数据假频。

图6展示的是对该数据一条inline线共偏移距射线束偏移成像结果，偏移输出网格深度采样间隔为2m，总深度6km，crossline采样间隔12.5m,输出范围为13.75km。图6(a)偏移基于无反假频波束形成，图6(b)和图6(c)基于不同权重的反假频波束形成，分别为wx＝wy＝1.0(见图6(b))和wx＝wy＝1.5(见图6(c))。无假频成像结果中整个剖面上随机噪音严重，尤其是浅层薄层结构被掩盖难于辨识。反假频后成像结果(见图6(b)和图6(c))随机噪音被明显压制，浅层水平结构辨识度增高。图7为选择该inline线中部分区域并用波形图放大显示其结构细节。在crossline1.8km、深度1.5km左右的断层(见图7(a)、7(c)和7(e))在反假频结果中可被明显辨识(见图7(b))。对比图6(b)和图6(c)，图6(c)中浅部噪音被进一步压制，浅部薄层结构分辨率进一步增强，但一些陡倾角结构明显被削弱(见图7(b)、(d)和(f))。

图8所示为对应图6成像结果的共成像点偏移距道集，每个成像点输出20个偏移距，最小偏移距200m，偏移距间隔200m。无假频道集中(见图8(a))随机噪音严重且分布于所有偏移距上，尤其是浅部和深部的弱反射信号难于辨识，使avo分析或基于偏移距道集的速度分析难于进行。经反假频后道集反射信号连续性和相干性明显增强(见图8(b)、(c)),反假频后小偏移距分辨率变化较小，大偏移距能量和分辨率均有所降低。

综合考虑偏移效果和效率，将反假频系数为(wx＝1.0，wy＝1.0)应用于体偏移，该套参数压制了部分偏移噪音，而且对陡倾角成像分辨率和能量保持较好，剩余噪音可通过叠后剖面去噪方法进一步去除。图9所示为该海洋数据体偏移结果，以图6中所示inline为中心的101条线输出，在测试线两侧各50条，总长度2.5km。所选深度分别为500m(图9(a)),3km(图9(b))和5km(图9(c))，可以看到浅层仍保留了一定偏移噪音，在3km和5km深度剖面上信噪比较高，偏移噪音得到有效压制。

本发明的三角反假频滤波算子直接在时域实时完成，在对地震数据进行因果和反因果积分后，可以直接与时域倾斜叠加波束形成结合完成，避免了将数据变换至频域进行滤波操作。同时，本发明中提出在反假频公式中加入权重系数，以对反假频程度进行控制。本发明提出的反假频方法是在叠前数据域完成，根据地表速度决定的慢度分量在数据域滤波，反假频后数据可进行压缩存储后重复使用，不受速度模型更新影响。