本发明是关于一种时移地震虚拟震源双向波场重构方法和系统，属于地震数据波处理技术领域。
背景技术：
时移地震技术是利用不同时间观测的地震有效信息的差异，进行储层监测、完善油藏管理方案、提高油气采收率的技术。时移地震的可重复性是长期油藏监测最重要的影响因素。时移地震的可重复性与采集中的仪器因素，环境、背景噪声等采集因素密切相关。其中仪器因素较易控制，而环境、背景噪声往往很难控制，是造成监测结果不可重复的主要原因。环境、背景噪声主要源于由震源、检波器和近地表条件的变化导致的时间延迟和相位、振幅的畸变。现有技术中，通常通过将震源和检波器埋入地下的方式来提升时移地震的可重复性。但对于沙漠环境，近地表发育溶岩结构的地区的成像结果尤为复杂，会引入多次波、面波、散射噪声，以及不同波型互相干扰，使上下行波场发生畸变。
为了解决上述问题，引入虚拟震源法(vs)对时移地震进行检测，虚拟震源法不依赖于近地表速度模型，其将上下行波场进行互相关重构，将地面震源重构到地下观测系统的位置，压制来自近地表仪器耦合和观测系统变化形成的噪声，在储层监测和时移地震处理中，提高了可重复性。但是，所有基于干涉的波场重构法仅涉及单向重构，即将震源或检波器一端重构到地下观测系统的位置，其虽然可以避免地表因素对波场的影响，但由于震源或检测器位于地下观测系统，而与其相对的检测器或震源则在地表，所以重构波场仍然需要经过地层结构。并且，为了得到高质量的成像效果，需要高密度观测系统，成本高昂。
技术实现要素：
针对上述现有技术的不足，本发明的目的是提供了一种时移地震虚拟震源双向波场重构方法和系统，其将震源和检波器都重构到地表以下的观测系统，并对重构数据插值化，降低了地层结构对波场的影响，以较低的成本实现高质量地震波场成像。
为实现上述目的，本发明提供了一种时移地震虚拟震源双向波场重构方法，包括以下步骤：s1.将预先获取的地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；s2.将经过重构的地震数据信号变换回时域；s3.对时域的地震数据信号进行共中心点叠加，得到虚拟震源激发波场
进一步，地震数据信号分为地下观测数据和地表观测数据，地下观测数据包括震源端的直达波和接收端的直达波；地表观测数据包括去除面波的反射波，将震源端的直达波、接收端的直达波和反射波转换至频率域。
进一步，步骤s2中震源端的波场重构公式为：
其中，x、y和b分别为震源、地表接收器和地下接收器的空间坐标；v(b|y；ω)是地下接收器b作为垂向虚拟震源时，地表接收器y接收到的波场；v*(b|y；ω)是v(b|y；ω)的复共轭；d*(x|b；ω)是地下接收器b接收到震源激发的直达波场；u(x|y；ω)是地表检波器y接收到的来自震源x的反射上行波；nx是积分面的法向量。
进一步，步骤s2中接收端的波场重构公式为：
其中，v(a|b；ω)是将地下接收器a视为虚拟源时，另一地下接收器b处的重构波场；k是地表震源发出的波数；d*(x|a)是地下检波器a接收到的地表震源x激发的波场；d*(b|x；ω)是将地下接收器b视为虚拟源时，震源x处接收到的重构波场。
进一步，对d*(x|a)和d*(b|x；ω)进行零相位白化滤波，经过滤波后的重构波场为：
其中，h为零相位白化滤波函数，*为复共轭。
进一步，步骤s3中将经过重构的地震数据信号变换回时域的公式为：
vh(a|b；t)＝fω→t{vh(a|b；ω)}
其中，fω→t表示频率域到时域的傅里叶变换。
进一步，地下观测数据的接收器少于地表观测数据的接收器，为了补偿地下的接收器覆盖次数不足，步骤s3中将地震数据信号转换为时域后，需对地震数据信号进行插值处理。
进一步，采用平面波分解插值方法对地震数据信号进行处理，插值写作最小二乘方程：
其中，vn是第n条道的经过重构的地震数据信号，pn,n-1是第n道与第n-2道之间的倾角域平面波分解算子，通过公式：
进一步，平面波分解插值方法，插值后得到的重构波场为：
其中，
本发明还公开了一种时移地震虚拟震源双向波场重构系统，包括：双向波场重构模块，用于将预先获取的地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；时域转换模块，用于将经过重构的地震数据信号变换回时域；共中心点叠加模块，用于对时域的地震数据信号进行共中心点叠加，得到虚拟震源激发波场。
本发明由于采取以上技术方案，其具有以下优点：
1、本发明使用稀疏观测系统，将震源和检波器都重构到地表以下的观测系统，并对重构数据插值化，降低了地层结构对波场的影响，以较低的成本实现高质量地震波场成像。
2、本发明将虚拟源波场重构方法从单向拓展到双向，进一步压制了复杂近地表干扰，提高了成像质量与可重复性；
3、本发明将倾角引导的地震道插值方法应用于波长重构，补偿了稀疏观测系统覆盖次数的不足，降低了成本。
附图说明
图1是本发明一实施例中时移地震虚拟震源双向波场重构方法的流程图；
图2是在复杂近地表背景下，时移地震虚波场检测的原理示意图；
图3是地表接收器y接收的地震波信号；
图4是地下接收器a和b接收的地震波信号；
图5是双向波场重构中的震源端的波场重构的原理示意图；
图6是双向波场重构中的接收端的波场重构的原理示意图；
图7是重构的接收波场中的下行直达波d(y|a)的示意图；
图8是震源x到地下检波点a的直达波d(x|a)按行偏移量排列的单道信号记录；
图9中图9(a)是经过双向波场重构后获得的单炮信号记录；图9(b)是使用平面波分解滤波器进行初始倾角估算后得到的单炮信号记录；图9(c)是经过插值后的单炮信号记录；
图10是经过双向波场重构和插值后的波场检测的原理示意图；
图11是地震波测试数据观测系统的布局图；
图12是本发明一实施例中采集到的三维数据体示意图；
图13中图13(a)是地表观测系统共中心点叠加剖面图，图13(b)是地下观测系统共中心点叠加剖面图，图13(c)是经过双向重构波场共中心点叠加剖面图；
图14为第一组勘探(第1-6组)和第二组勘探(第7-13组)的叠加剖面图对比，左侧是第一组勘探剖面图，右侧是第二组勘探剖面图，图14(a)是地表观测系统共中心点叠加剖面图，图14(b)是地下观测系统共中心点叠加剖面图，图14(c)是经过双向重构波场共中心点叠加剖面图；
图15是地表观测系统、地下观测系统和经过双向重构波场系统的均方根均值的示意图。
具体实施方式
为了使本领域技术人员更好的理解本发明的技术方向，通过具体实施例对本发明进行详细的描绘。然而应当理解，具体实施方式的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。在本发明的描述中，需要理解的是，所用到的术语仅仅是用于描述的目的，而不能理解为指示或暗示相对重要性。
实施例一
本实施例公开了一种时移地震虚拟震源双向波场重构方法，如图1所示，包括以下步骤：
s1.获取地震数据信号；
s2.将地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；
s3.将经过重构的地震数据信号变换回时域；
s4.对地震数据信号进行共中心点叠加，得到虚拟震源激发波场。
本实施例中方法将震源和检波器都重构到地表以下的观测系统，并对重构数据插值化，降低了地层结构对波场的影响，以较低的成本实现高质量地震波场成像。
地震数据信号分为地下观测数据和地表观测数据，地下观测数据包括震源端的直达波和接收端的直达波；地表观测数据包括去除面波的反射波，将震源端的直达波、接收端的直达波和反射波转换至频率域。
反射波地震干涉法是不同接收器接收到同一震源发出的地震信号，并将各个地震信号进行互相关，然后对每个相关震源进行积分得到的地震波场图。如图2所示，图2为在复杂近地表背景下，时移地震虚波场检测的原理示意图，图中星型图案代表震源x，倒三角型图案代表地表接收器y。a和b都表示设置在地下的接收器。此处的接收器一般为检波器，用于接收地震波信号。图中最下端的圆型图案代表目标反射体。震源生成的下行直达波d(x|b)先被地下接收器b接收并记录，经过目标反射体反射后转换为上行波u(x|a)被地下接收器a接收并记录，上行波u(x|a)继续向地表传播，直到被地表接收器y接收并记录，生成原始的地表记录u(x|y)。图3是地表接收器y接收的地震波信号。图3中的横轴为炮号，单炮记录中的每一道按照炮检距排列，纵轴为时间。图4是地下接收器a和地下接收器b接收的地震波信号。图4与图3类似，图4中的横轴为炮号，单炮记录中的每一道按照炮检距排列，纵轴为时间。将图4中的地震信号图与图3中的地震信号图进行比较后可以看出，地下接收器a和b记录的信号具有更高的信噪比和更好的反射效果。其中，原始的近偏移距波场被视为直达下行波场(偏移距<30m，时间<100ms)，因为早期的波场主要由直达p波组成。类似的，时间上较晚到达的波场(时间>300ms)被视为上行反射波场。基于图3和图4的分析，需要将位于地表的震源和接收器的波场重构至地表以下，以避免地表的情况和地层结构对地震波场的影响。
本实施例中步骤s2中震源端的波场重构的原理如图5所示，地表震源发出的，地表接收器y接收到的地震波信号数据u(x|y)与地下接收器b接收到的下行直达波d(x|b)互相关，并进一步将多个单炮记录叠加，以得到重构震源波场v(b|y)。
v(b|y)满足以下公式：
其中，ω是角频率，x、y和b分别为震源、地表接收器和地下接收器的空间坐标；v(b|y；ω)是地下接收器b作为垂向虚拟震源时，地表接收器y接收到的波场；v*(b|y；ω)是v(b|y；ω)的复共轭；d*(x|b；ω)是地下接收器b接收到震源激发的直达波场；u(x|y；ω)是地表检波器y接收到的来自震源x的反射上行波；nx是积分面的法向量。
公式(1)在远场地震的情况下，可以近似的简化为如下形式：
v(b|y；ω)-v*(b|y；ω)≈2ik∫surfaced*(x|b；ω)u(x|y；ω)d2x(2)
k是地表震源发出的波数。远场地震为震中距大于1000公里的地震。尽管震源波场重构将地下接收器重构为虚拟源v(b|y；ω)，但仍然无法完全摆脱复杂近地表的干扰，因为地震信号到达地表检波器y时仍不可避免地穿过近地表层。
故需要将接收端的波场也需要重构为地下接收器波场，如图6所示，图6为接收端的波场重构的原理示意图。将以地下接收器b为虚拟震源的信号v(b|y)转换为以地表接收器y为虚拟震源的信号v(y|b)，接着与地下检波器a接受到的来自虚拟震源y的直达波d(y|a)互相关，并将多个以y为震源的单炮信号叠加，生成重构波场v(a|b)。
接收端的波场重构v(a|b；ω)满足以下公式：
v(a|b；ω)-v*(a|b；ω)≈2ik∫buriedd*(y|a；ω)u(y|b；ω)d2y(3)
其中，v(a|b；ω)是将地下接收器a视为虚拟源时，另一地下接收器b处的重构波场；v*(a|b；ω)是v(a|b；ω)的复共轭；k是地表震源发出的波数；d*(y|a)是地下检波器a接收到的虚拟震源y激发的波场；u(y|b；ω)是地下接收器b接收到的来自虚拟震源y的反射上行波。；d*(y|a)是地下检波器a接收到的虚拟震源y激发的波场，但y本身无法激发地震波，所以需要将真实震源x与地表接收器y设置的尽可能接近，如图7所示。图7是重构的接收波场中的下行直达波d(y|a)的示意图。真实震源x与地表接收器y非常接近，因此d(x|a)≈d(y|a)。图8是震源x到地下检波点a的直达波d(x|a)按行偏移量排列的单道信号记录。从图8中可以看出直达波的波形非常复杂，体现了近地表结构的复杂性。图8中由于地下检波点a是单分量检波器，且检测地震波的垂直分量，故小角度的波呈现出大幅度的振幅变化，而大角度的波则相对较小的振幅变化，因为这是垂直分量记录。
基于图7中的近似，将公式(2)带入公式(3)得出：
将d(x|a；ω)带入u(y|x；ω)，公式(4)可以简化为：
v(a|b；ω)≈-4k2∫surface∫buriedd*(x|a)u(y|x；ω)d*(b|x；ω)d2xd2y(5)
其中，v(a|b；ω)是将地下接收器a视为虚拟源时，另一地下接收器b处的重构波场；k是地表震源发出的波数；d*(x|a)是地下检波器a接收到的地表震源x激发的波场；d*(b|x；ω)是将地下接收器b视为虚拟源时，震源x处接收到的重构波场。在共道集中，可以将地表数据与地下检波器记录的直达波互相关，实现双向波场重构。与经典的地震波场重构方法不同，公式(5)完全由数据驱动，受地表层以及人工因素的影响较小。图9中图9(a)是经过双向波场重构后获得的单炮信号记录，与图3和图4中的单炮信号记录相比，图9(a)中的单炮信号记录的信噪比更高，尤其是对于浅层目标反射效果更好，从而减少了近地表复杂性对地震波场的影响。
在频率域中，与原始地表数据相比，双向重构波场v(a|b)的频谱不够精确。为了得到更加精确的对比双向重构波场与非重构波场，需要对d(x|a)和d(b|x)应用零相位白化滤波函数h，以保证双向重构波场的频谱与u(y|x)一致。将零相位白化滤波函数h加入公式(5)可得到最终重构结果：
vh(a|b；ω)≈-4k2∫surface∫buried[hd*(x|a；ω)]*u(y|x；ω)[hd*(b|x；ω)]*d2xd2y(6)
其中，h为零相位白化滤波函数，*为复共轭。
地下观测数据的接收器少于地表观测数据的接收器，为了补偿地下的接收器覆盖次数不足，步骤s3中将地震数据信号转换为时域后，需对地震数据信号进行插值处理。
vh(a|b；ω)包含横向变化较小的水平反射层，适合使用基于平面波分解滤波器的倾角引导的地震道插值方法。
平面波分解滤波器(pwd)是一种评估地震记录倾角的方法，它通过将上一道沿主同相轴相移，保幅地预测下一道。然后，通过最小化预测误差来计算估计的主坡度。通过形状正则化保证了重构波场中估计的倾角变化平稳。由于pwd通常应用于时域，故需要将经过下面的公式将重构的地震数据信号变换回时域：
vh(a|b；t)＝fω→t{vh(a|b；ω)}(7)
其中，fω→t表示频率域到时域的傅里叶变换。
采用平面波分解插值方法对地震数据信号进行处理，插值写作最小二乘方程：vhp≈δ，或者
其中，vn是第n条道的经过重构的地震数据信号，pn,n-1是第n道与第n-2道之间的倾角域平面波分解算子，通过公式：
公式：
由于对角矩阵的逆非常容易求解，因此可以直接求解方程(9)来获得倾角场。图9(b)是使用平面波分解滤波器进行初始倾角估算后得到的单炮信号记录；对倾角场进行了平滑处理以形成扩展的倾角场
另一种重构
使用结构约束插值，插值后的重构波场为：
其中，
图10是经过双向波场重构和插值后的波场检测的原理示意图。通过该原理示意图，可以看出双向波场干涉互相关将震源和观测系统均“布设”于地下，避免了复杂近地表干扰。
实施例二
本实施例用沙漠环境中的时移地震数据，对实施例一中的方法进行了进一步的说明和测试。
在沙漠地区，19个月内进行了13次重复二维勘探，前6次(s1–s6)在3个月完成，中断了17个月后又在一周内完成后7次勘探(s7–s13)。13次勘探均使用mertz26可控震源，震源重复激发的精度达1m。地下检波器的布局方式如图11所示，图中三角型图案为80个地下检波器，点状图案为地表的震源。所有的震源组成9条炮线。中间炮线附近还布设有300个地表检波器，包含300个地表检波器的测线紧邻第五条炮线，即9条炮线的中间的一条炮线，道间距为7.5m，最小偏移距小于5米。近地表为厚度从几米到几十米不等的砂层，近地表以下是倾角小于5°的层状结构。数据的成像质量和重复性受到近地表结构干扰。
本实施例中的步骤一获取地震数据信号分三步处理野外数据：第一步包括九条炮线多模式正交求和，以去除过强的爆破信号、随机能量和散射噪声。利用高密度炮点，将共炮道级整合为共检波点道级，用fk滤波压制面波噪声。第二步是对数据进行双向波场重构。用自动算法拾取60ms窗口内的直达波，并根据直达波求出重构算子d(x|a)和d(b|x)。第三步，对数据进行静校正和动校正以及振幅补偿、顶切和共中心点叠加。振幅补偿是对每道的5个500ms窗口逐个补偿。共中心点叠加包括算数求和以及对不同覆盖次数的道集进行能量平均。用于对比的两种非重构的传统共中心点叠加方法只利用了地表波场u(y|x)和地下波场v(y|b)，未经过第二步重构处理。
图12是本实施例中采集到的三维数据体示意图，其中，x和y轴分别表示地表炮点(2700)和检波点(300/80)，xz平面表示共道集，yz平面表示共炮集。本实施例中使用矩阵转置的方法进行单炮记录与单道记录之间的转换。
图13中图13(a)是地表观测系统共中心点叠加剖面图，图13(b)是地下观测系统共中心点叠加剖面图，图中箭头标出的层为目标反射层。图13(a)中，地表测线叠加剖面图上的信噪比最低，尤其剖面图的右侧受到地表沙漠环境的影响，成像效果非常差。如图13(b)所示，图13(b)的地下观测系统共中心点叠加剖面图，即使不重构，信噪比和连续性都比图13(a)中要好，这是因为地下检波器部分避免了地表层对地震波的影响。图13(c)是经过双向波场重构后的叠加剖面图，图13(c)中剖面图的成像质量相比图13(a)的图13(b)明显提升，即使是较浅的目标层也完全没有受到地表砂层的干扰，波场的连续性增强。这些结果证明了本实施例中方法的鲁棒性，并且这些成像结果还可以进一步应用于地层速度分析和偏移。
除了图像质量，可重复性也是衡量地震监测质量的一个重要因素。为对比本实施例中测试系统的可重复性，从间隔一年的两组勘探(地下测线不变)中随机抽取一道。如图14所示，图14为第一组勘探(第1-6组)和第二组勘探(第7-13组)的叠加剖面图对比，左侧是第一组勘探剖面图，右侧是第二组勘探剖面图，图14(a)是地表观测系统共中心点叠加剖面图，图14(b)是地下观测系统共中心点叠加剖面图，图14(c)是经过双向重构波场共中心点叠加剖面图。从图14(a)中可以看出，剖面图左面部分和右面部分存在很明显的差异，也就是说两组勘探的结果重复性很差，尤其是图中箭头标记出的目标层。图14(b)的剖面图左面部分和右面部分也存在差异，但相对于图14(a)，图14(b)的剖面图中波场的连续性更好，尤其是箭头标记处的目标层。因为地下测线在接收端减少了部分近地表干扰。但由于震源端仍经过近地表，所以仍存在部分不一致。图14(c)的剖面图左面部分和右面部分的差异很小，可见本实施例中，将地表震源和检波器都重基准到近地表以下，可显著提高检测的重复性。
图15是地表观测系统、地下观测系统和经过双向重构波场系统的均方根均值的示意图。可以通过在目标层窗口内计算归一化均方根(nrms)的方法来定量评价测试的重复性。nrms是给定窗口内两个地震剖面的均方根占均值的百分比。nmrs对地震数据的细微变化十分敏感。分别对13次勘探结果计算nrms，每个共中心点的道集可求得78个nrms值。三种方法的nrms直方图如图15所示。地表观测系统为双峰分布，峰值分别出现在35％和100％。两组勘探间的nrms值较大，两个分离的峰顶具有不同的形状，这表明地表观测系统重复性问题差。对于地下观测系统同样呈双峰分布，峰值出现在约20％和70％处。相对于地表观测系统分布更加集中，但两组勘探间仍存在不一致的问题。经过双向重构波场系统重构数据nrms呈现单峰分布，峰值位于32％。可见，相比于传统虚拟源法，经过双向重构波场系统显著提升了成像质量和可重复性。
实施例三
基于相同的发明构思，本实施例公开了一种时移地震虚拟震源双向波场重构系统，包括：
数据获取模块，用于获取地震数据信号；
双向波场重构模块，用于将地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；
时域转换模块，用于将经过重构的地震数据信号变换回时域；
共中心点叠加模块，用于对地震数据信号进行共中心点叠加，得到虚拟震源激发波场。
上述内容仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围为准。
1.一种时移地震虚拟震源双向波场重构方法，其特征在于，包括以下步骤：
s1.将预先获取的地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；
s2.将经过重构的所述地震数据信号变换回时域；
s3.对时域的地震数据信号进行共中心点叠加，得到虚拟震源激发波场。
2.如权利要求1所述的时移地震虚拟震源双向波场重构方法，其特征在于，所述地震数据信号分为地下观测数据和地表观测数据，所述地下观测数据包括震源端的直达波和接收端的直达波；所述地表观测数据包括去除面波的反射波，将震源端的直达波、接收端的直达波和反射波转换至频率域。
3.如权利要求2所述的时移地震虚拟震源双向波场重构方法，其特征在于，所述步骤s2中震源端的波场重构公式为：
其中，x、y和b分别为震源、地表接收器和地下接收器的空间坐标；v(b|y；ω)是地下接收器b作为垂向虚拟震源时，地表接收器y接收到的波场；v*(b|y；ω)是v(b|y；ω)的复共轭；d*(x|b；ω)是地下接收器b接收到震源激发的直达波场；u(x|y；ω)是地表检波器y接收到的来自震源x的反射上行波；nx是积分面的法向量。
4.如权利要求3所述的时移地震虚拟震源双向波场重构方法，其特征在于，所述步骤s2中接收端的波场重构公式为：
其中，v(a|b；ω)是将地下接收器a视为虚拟源时，另一地下接收器b处的重构波场；k是地表震源发出的波数；d*(x|a)是地下检波器a接收到的地表震源x激发的波场；d*(b|x；ω)是将地下接收器b视为虚拟源时，震源x处接收到的重构波场。
5.如权利要求4所述的时移地震虚拟震源双向波场重构方法，其特征在于，对所述d*(x|a)和d*(b|x；ω)进行零相位白化滤波，经过滤波后的重构波场为：
其中，h为零相位白化滤波函数，*为复共轭。
6.如权利要求1-5任一项所述的时移地震虚拟震源双向波场重构方法，其特征在于，步骤s3中将经过重构的所述地震数据信号变换回时域的公式为：
vh(a|b；t)＝fω→t{vh(a|b；ω)}
其中，fω→t表示频率域到时域的傅里叶变换。
7.如权利要求2-5任一项所述的时移地震虚拟震源双向波场重构方法，其特征在于，所述地下观测数据的接收器少于地表观测数据的接收器，为了补偿地下的接收器覆盖次数不足，所述步骤s3中将所述地震数据信号转换为时域后，需对所述地震数据信号进行插值处理。
8.如权利要求7所述的时移地震虚拟震源双向波场重构方法，其特征在于，采用平面波分解插值方法对所述地震数据信号进行处理，所述插值写作最小二乘方程：
其中，vn是第n条道的经过重构的地震数据信号，pn,n-1是第n道与第n-2道之间的倾角域平面波分解算子，通过公式：
9.如权利要求8所述的时移地震虚拟震源双向波场重构方法，其特征在于，所述平面波分解插值方法，插值后得到的重构波场为：
其中，
10.一种时移地震虚拟震源双向波场重构系统，其特征在于，包括：
双向波场重构模块，用于将预先获取的地震数据信号转换至频率域，对震源端和接收端的波场进行重构，将震源和接收器都重构到地表以下；
时域转换模块，用于将经过重构的所述地震数据信号变换回时域；
共中心点叠加模块，用于对时域的地震数据信号进行共中心点叠加，得到虚拟震源激发波场。
技术总结