本发明是关于一种基于小波互相关时移地震虚拟震源波场重构方法和系统，属于地震数据波处理技术领域。
背景技术：
时移地震技术是利用不同时间观测的地震有效信息的差异，进行储层监测、完善油藏管理方案、提高油气采收率的技术。时移地震的可重复性是长期油藏监测最重要的影响因素。时移地震的可重复性与采集中的仪器因素，环境、背景噪声等采集因素密切相关。其中仪器因素较易控制，而环境、背景噪声往往很难控制，是造成监测结果不可重复的主要原因。环境、背景噪声主要源于由震源、检波器和近地表条件的变化导致的时间延迟和相位、振幅的畸变。现有技术中，通常通过将震源和检波器埋入地下的方式来提升时移地震的可重复性。但对于沙漠环境，近地表发育溶岩结构的地区的成像结果尤为复杂，会引入多次波、面波、散射噪声，以及不同波型互相干扰，使上下行波场发生畸变。
为了解决上述问题，引入虚拟震源法(vs)对时移地震进行检测，虚拟震源法不依赖于近地表速度模型，其将上下行波场进行互相关重构，将地面震源重构到地下观测系统的位置，压制来自近地表仪器耦合和观测系统变化形成的噪声，在储层监测和时移地震处理中，提高了可重复性。然而，由于有限的震源孔径和野外预处理的结果不能完全满足虚拟震源的假设条件，因此仍然存在重复性差和图像质量差的问题。尤其，在强烈非均质的上覆层中，下行波场往往被近地表反射和自由表面多次波的上行波干扰，使二者难以分离。一旦下行的p波受到污染，就无法有效抑制包含多次波、散射和面波的噪声和串扰，导致传统虚拟震源叠加的信噪比低于实际震源叠加，使成像质量下降。
技术实现要素：
针对上述现有技术的不足，本发明的目的是提供了一种基于小波互相关时移地震虚拟震源波场重构方法和系统，在传统虚拟震源方法处理流程中，加入小波域虚拟震源波场重构，对小波系数进行时间尺度域ts或时间尺度波数域tsk滤波，从而有效压制噪声，并对非平稳信号实现高分辨率分离。该方法无需近地表模型即可有效地抑制伪像，优化叠加成像的结果。
为实现上述目的，本发明提供了一种基于小波互相关时移地震虚拟震源波场重构方法，包括以下步骤：s1.通过若干检波器获得垂直地震剖面vsp地震数据，将vsp地震数据进行上行波和下行波分离；s2.将分离后的上行波和下行波都进行小波变换，选择其中一个检波器作为虚拟震源，把检波器的下行波小波系数与其它检波器的上行波小波系数互相关，得到小波互相关系数；s3.对小波互相关系数进行时间尺度域ts或时间尺度波数域tsk滤波；s4.将经过滤波的小波互相关系数逆变换回时空域，得到虚拟震源道集；s5.将虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
进一步，步骤s1中垂直地震剖面vsp地震数据采用下式表示：
其中，
进一步，步骤s1中垂直地震剖面vsp地震数据分别通过去噪和加窗获得上行波和下行波，上行波和下行波的公式为：
u(rb|rs；t)＝up(rb|rs；t) um(rb|rs；t) ush(rb|rs；t)...
d(rb|rs；t)＝dp(rb|rs；t) dm(rb|rs；t) dsh(rb|rs；t)...
其中，d(rb|rs；t)是检波器rb接收到的下行波场，up(rb|rs；t)、dp(rb|rs；t)分别是直达p波的上行波场和下行波场，um(rb|rs；t)、dm(rb|rs；t)分别为多次波近地表反射的上行波场和下行波场，ush(rb|rs；t)、dsh(rb|rs；t)分别为s波的上行波场和下行波场。
进一步，步骤s2中进行小波变换为：对于一个给定的地震信号d(t)采用以下公式进行小波变换：
进一步，步骤s2中小波互相关系数采用计算公式为：
其中，wv(α，τ)为小波互相关系数，
进一步，小波互相关系数满足如下条件：
wv(α，ω)＝αv(ω)|φ(αω)|2
|wv(α，ω)|＝α|v(ω)||φ(αω)|2
∠wv(α，ω)＝∠v(ω)
wv(α，ω)是小波尺度和频域中的小波互相关系数，v(ω)是信号互相关谱，φ(αω)是母小波的频谱，∠表示相位谱，ω是角频率，α为小波尺度。
进一步，步骤s3中时间尺度空间域的三维滤波方程为：
wvo(α，τ，x)＝h(α，τ，x)*α，τwvi(α，τ，x)
其中，wv0(α，τ，x)是时间尺度空间域中小波互相关经过滤波的小波互相关系数；h(α，τ，x)是时间尺度空间域中的滤波函数；α，τwvi(α，τ，x)是时间尺度空间域中的小波互相关原始信号。
进一步，对步骤s3中时间尺度空间域的三维滤波方程进行二维傅立叶变换，得到时间频率波数域，所述时间频率波数域的三维滤波方程为：
其中，wv0(α，ω，k)是滤波结果在时间频率波数域中小波互相关经过滤波的小波互相关系数，vv(ω，k)是频率波数域中虚拟震源道集的互相关谱，p(α，ω，k)是依赖于尺度因子α的滤波函数；
进一步，小波互相关系数的逆变换公式为：
其中，d(t)为经过处理后的地震信号，
基于相同的发明构思，本发明还公开了一种基于小波互相关时移地震虚拟震源波场重构系统，包括：上下行波分离模块，用于将垂直地震剖面vsp地震数据进行上行波和下行波分离；小波互相关系数获取模块，用于将上行波和下行波都进行小波变换，选择一检波器作为虚拟震源，把检波器的下行波小波系数，与其它检波器的上行波小波系数互相关，得到小波互相关系数；滤波模块，用于对小波互相关系数进行时间尺度域、时间尺度波数域或时间空间域滤波；逆变换模块，用于将经过滤波的小波互相关系数逆变换回时空域，得到虚拟震源道集；共中心点叠加模块，用于将虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
本发明由于采取以上技术方案，其具有以下优点：
1、在传统虚拟震源方法处理流程中，加入小波域虚拟震源波场重构，对小波系数进行时间尺度域ts或时间尺度波数域tsk滤波，从而有效压制噪声，并对非平稳信号实现高分辨率分离，该方法无需近地表模型即可有效地抑制伪像，优化叠加成像的结果。
2、本发明将地震干涉方法从一维时间域拓展到了多维ts域和tsk域，提高信噪比，使依赖于频率和波数的噪声和信号能够更好地分离。
3、相比于常规虚拟震源方法，本发明能够有效抑制与近地表结构相关的噪声，大幅提升时移地震的可重复性。
附图说明
图1是本发明一实施例中基于小波互相关时移地震虚拟震源波场重构方法的流程图；
图2是经过共道集噪声处理的单炮信号图，其中，上方的虚线上方为p波的下行波场，下方的虚线下方是p波的上行波场；
图3中图3(a)是某一虚拟单炮信号的下行波场；图3(b)是该下行波场对应的频率-波数谱；
图4中图4(a)是经过ts滤波后的fk域的下行波场；图4(b)是过ts滤波后的下行波场对应的频率-波数谱，图4中箭头标出的部分与图3中相应位置的图像相比，图4(a)的图像清晰度提高，图4(b)的波形噪音较低；
图5中图5(a)是仅经过ts滤波的虚拟单炮信号图；图5(b)是经过ts滤波和tsx三维滤波的虚拟单炮信号图；
图6是对tsx三维滤波进行二维傅立叶变换后得到的tsk三维数据体；
图7是进行sk滤波后得到的三维数据体；
图8中图8(a)是采用传统虚拟震源法经过共中心点叠加处理的波场剖面图，图8(b)是采用本发明一实施例中方法经过共中心点叠加处理的波场剖面，将图8(a)和图8(b)中箭头标出的部分进行比较，可见图8(b)中图像更加清晰，噪音较低；
图9是13次时移地震的二维剖面形成的三维数据体，图9(a)非虚拟震源法得到的三围数据体；图9(b)现有技术中虚拟震源法得到的三围数据体；图9(c)本实施例中虚拟震源法得到的三围数据体。
具体实施方式
为了使本领域技术人员更好的理解本发明的技术方向，通过具体实施例对本发明进行详细的描绘。然而应当理解，具体实施方式的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。在本发明的描述中，需要理解的是，所用到的术语仅仅是用于描述的目的，而不能理解为指示或暗示相对重要性。
实施例一
本实施例公开了一种基于小波互相关时移地震虚拟震源波场重构方法，如图1所示，包括以下步骤：
s1.通过若干检波器获得垂直地震剖面vsp地震数据，并将垂直地震剖面vsp地震数据进行上行波和下行波分离。
其中，垂直地震剖面vsp地震数据采用下式表示：
其中，
垂直地震剖面vsp地震数据分别通过去噪和加窗获得上行波和下行波，上行波和下行波的公式为：
u(rb|rs；t)＝up(rb|rs；t) um(rb|rs；t) ush(rb|rs；t)...
d(rb|rs；t)＝dp(rb|rs；t) dm(rb|rs；t) dsh(rb|rs；t)...
其中，d(rb|rs；t)是检波器rb接收到的下行波场，up(rb|rs；t)、dp(rb|rs；t)分别是直达p波的上行波场和下行波场，um(rb|rs；t)、dm(rb|rs；t)分别为多次波近地表反射的上行波场和下行波场，ush(rb|rs；t)、dsh(rb|rs；t)分别为s波的上行波场和下行波场。
s2.将上行波和下行波都进行小波变换，选择一检波器作为虚拟震源，把检波器的下行波小波系数，与其它检波器的上行波小波系数互相关，得到小波互相关系数。
其中，进行小波变换步骤具体为：对于一个给定的地震信号d(t)采用以下公式进行小波变换：
其中，
通过上行波和下行波的小波系数互相关计算小波互相关系数，小波互相关系数采用如下公式计算：
其中，wv(α，τ)为小波互相关系数，
为了保证小波尺度和频域中的小波的同相，小波互相关系数满足如下条件：
wv(α，ω)＝αv(ω)|φ(αω)|2
|wv(α，ω)|＝α|v(ω)||φ(αω)|2
∠wv(α，ω)＝∠v(ω)
wv(α，ω)是小波尺度和频域中的小波互相关系数，v(ω)是信号互相关谱，φ(αω)是母小波的频谱，∠表示相位谱，ω是角频率，α为小波尺度。该方法仅对振幅谱滤波，保留原始相位，并从虚拟震源中提取出格林函数的运动学特征，是一种保相技术。
s3.对小波互相关系数进行时间尺度域ts或时间尺度波数域tsk滤波。
其中，时间尺度空间域的三维滤波方程为：
wvo(α，τ，x)＝h(α，τ，x)*α，τwvi(α，τ，x)
其中，wv0(α，τ，x)是时间尺度空间域中小波互相关经过滤波的小波互相关系数；h(α，τ，x)是时间尺度空间域中的滤波函数；α，τwvi(α，τ，x)是时间尺度空间域中的小波互相关原始信号。
其中，滤波函数h(α，τ，x)旨在改善同相轴识别、处理非平稳特性并压制虚拟震源法处理过程中产生的相干噪声。三维滤波函数h(α，τ，x)表示为：
对时间尺度空间域的三维滤波方程进行二维傅立叶变换，得到时间频率波数域，所述时间频率波数域的三维滤波方程为：
其中，wv0(α，ω，k)是滤波结果在时间频率波数域中小波互相关经过滤波的小波互相关系数，v(ω，k)是频率波数域中虚拟震源道集的互相关谱，p(α，ω，k)是依赖于尺度因子α的滤波函数；
s4.将经过滤波的小波互相关系数逆变换回时空域，得到虚拟震源道集。
其中，小波互相关系数的逆变换公式为：
其中，d(t)为经过处理后的地震信号，
s5.虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
实施例二
本实施例以沙漠环境中的时移地震实际数据对实施例一中的方法进行了测试。具体过程如下：
在19个月内进行了13次重复二维勘探，前6次(s1-s6)在3个月完成，中断了17个月后又在一周内完成了后7次勘探(s7-s13)。13次勘探均使用了mertz26可控震源，震源重复激发的精度达1m。平均间隔30m的80个检波器和水听器固定在最深达50m的地下，组成一条二维测线。为了有效去除线性和散射噪声，实现最佳照明，进行了密集的三维面激发，将炮距设置为7.5×7.5m。偏移距范围为0到2400m。目标层位于地下约2000米处，上覆有多套速度差异较大的地层。近地表为厚度从几米到几十米不等的砂层，在近地表以下是简单的层状结构，反射界面的倾角均小于5°。
地震数据预处理时，先自动拾取初至的波形，再通过初至的波形的位置截取直达波的时窗用于下行波场近似。图2表示经过共道集噪声处理的单炮信号，上方虚线上面的数据无需处理，可直接用于p波的下行波场近似，对应30m或100ms以内的小偏移距。下方虚线下面的数据反射到达较晚，大于300ms，对应大偏移距，用于p波的上行波场近似。由于该区域地面情况复杂，下行波场往往包含来自近地表的干扰。
如图3所示，图3(a)是某一虚拟单炮信号的频率-波数域即fk域的下行波场；图3(b)是该下行波场对应的频率-波数谱。图3(a)中浅色区域为散射噪声，被浅色区域包裹的深色区域表示来自目标层的强反射。在fk域中的下行波场中的能量明显失真，因此，需进行适当滤波。与常规虚拟震源法相比，小波互相关和多通道组合滤波能将多维空间中的反射波信号从背景散射噪声中分离，提高信噪比，增加了反射波信号。本实施例中对每个炮集的上行波和下行波进行小波互相关。进行小波互相关时，由于在有效带宽内，目标储层的反射信号出现在约1.2-1.4s处。首先通过匹配目标储层的反射信号设置参数，本实施例中有效带宽优选为20-50hz，参数优选α＝20hz，τ＝0.2s和η＝5％的滑动二维窗口内的最大值。然后，将滤波函数与图3a中所示的小波系数振幅谱进行卷积。
经过小波互相关滤波的结果如图4所示，图4(a)是经过ts滤波后的fk域的下行波场；图4(b)是过ts滤波后的下行波场对应的频率-波数谱。在互相关计算之前，去除散射噪声能够进一步提升整体的信噪比。从图4(a)和图4(b)中可以看出经过滤波后，背景噪声明显减弱，但是浅层(＜0.6s)的高频成分(＞60hz)的信号被衰减，从而降低分辨率。
随后，采用多分量滤波器处理被面波污染的虚拟单炮信号。将经过ts滤波生成的虚拟单炮信号变换到时间尺度空间域上，进行tsx三维滤波，tsx三维滤波的结果如图5所示。图5(a)是仅经过ts滤波的虚拟单炮信号图；图5(b)是经过ts滤波和tsx三维滤波的虚拟单炮信号图。从图5(a)和图5(b)的对比中可以看出，经过tsx三维滤波的虚拟单炮信号强度明显提高，细节特征也更加明显，信噪比得到有效提高。
图6为对tsx三维滤波进行二维傅立叶变换后得到的tsk三维数据体。通过在tsk三维数据体中调整小波互相关系数可以获得特定的时间和频率分辨率。tsk滤波器包含时间信息，因而能够处理小波的非平稳性。当滤波器时间不变时，此方法简化为sk滤波，如图7所示。由此可见，即使使用简单的tsk滤波器，小波互相关也能够消除残留的面波，这是传统虚拟震源法所不能实现的。
图8中，图8(a)是采用传统虚拟震源法经过共中心点叠加处理的波场剖面图，图8(b)是采用本实施例中方法经过共中心点叠加处理的波场剖面。从图8(a)和图8(b)中可以看出，图8(b)的波场信号的信噪比和同相轴连续性的得到了显著提升。
图9是13次时移地震的二维剖面形成的三维数据体，图9(a)非虚拟震源法得到的三围数据体；图9(b)现有技术中虚拟震源法得到的三围数据体；图9(c)本实施例中虚拟震源法得到的三围数据体。从图9(a)、图9(b)和图9(c)中可以看出，经过本实施例中方法处理的时移地震波场具有更好的重复性。
本实施例中方法的优点在于能够有效消除强振幅噪声，特别是大偏移距数据中常出现的面波和陆上地震数据中的近地表溶岩产生的散射。不同的信号分解方法相互组合可以有效去除噪声。例如，小倾角的简单层状结构横向变化较小，ts滤波加tsx阈值是理想的选择。而对于包含断层或盐丘的地质结构，各向异性较强，存在潜在散射体，则需将信号转换到自适应tsk域或tk进行滤波。通用三维滤波比二维滤波更贵，但能更有效的解决高维度地震数据中的非平稳性。
实施例三
基于相同的发明构思，本实施例公开了一种基于小波互相关时移地震虚拟震源波场重构系统，包括：
上下行波分离模块，用于将垂直地震剖面vsp地震数据进行上行波和下行波分离；
小波互相关系数获取模块，用于将上行波和下行波都进行小波变换，选择一检波器作为虚拟震源，把检波器的下行波小波系数，与其它检波器的上行波小波系数互相关，得到小波互相关系数；
滤波模块，用于对小波互相关系数进行时间尺度域、时间尺度波数域或时间空间域滤波；
逆变换模块，用于将经过滤波的小波互相关系数逆变换回时空域，得到虚拟震源道集；共中心点叠加模块，用于将虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
上述内容仅为本申请的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应以权利要求的保护范围为准。
1.一种基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，包括以下步骤：
s1.通过若干检波器获得垂直地震剖面vsp地震数据，将所述vsp地震数据进行上行波和下行波分离；
s2.将分离后的上行波和下行波都进行小波变换，选择其中一个检波器作为虚拟震源，把所述检波器的下行波小波系数与其它检波器的上行波小波系数互相关，得到小波互相关系数；
s3.对小波互相关系数进行时间尺度域ts或时间尺度波数域tsk滤波；
s4.将经过滤波的小波互相关系数逆变换回时空域，得到虚拟震源道集；
s5.将所述虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
2.如权利要求1所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述步骤s1中垂直地震剖面vsp地震数据采用下式表示：
其中，
3.如权利要求2所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述步骤s1中垂直地震剖面vsp地震数据分别通过去噪和加窗获得上行波和下行波，上行波和下行波的公式为：
u(rb|rs；t)＝up(rb|rs；t) um(rb|rs；t) ush(rb|rs；t)...
d(rb|rs；t)＝dp(rb|rs；t) dm(rb|rs；t) dsh(rb|rs；t)...
其中，d(rb|rs；t)是检波器rb接收到的下行波场，up(rb|rs；t)、dp(rb|rs；t)分别是直达p波的上行波场和下行波场，um(rb|rs；t)、dm(rb|rs；t)分别为多次波近地表反射的上行波场和下行波场，ush(rb|rs；t)、dsh(rb|rs；t)分别为s波的上行波场和下行波场。
4.如权利要求1-3任一项所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述步骤s2中进行小波变换为：对于一个给定的地震信号d(t)采用以下公式进行小波变换：
其中，
5.如权利要求4所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，步骤s2中小波互相关系数计算公式为：
其中，wv(α，τ)为小波互相关系数，
6.如权利要求5所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述小波互相关系数满足如下条件：
wv(α，ω)＝αv(ω)|φ(αω)|2
|wv(α，ω)|＝α|v(ω)||φ(αω)|2
∠wv(α，ω)＝∠v(ω)
wv(α，ω)是小波尺度和频域中的小波互相关系数，v(ω)是信号互相关谱，φ(αω)是母小波的频谱，∠表示相位谱，ω是角频率，α为小波尺度。
7.如权利要求5或6所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述步骤s3中时间尺度空间域滤波采用三维滤波方程：
wvo(α，τ，x)＝h(α，τ，x)*α，τwvi(α，τ，x)
其中，wv0(α，τ，x)是时间尺度空间域中小波互相关经过滤波的小波互相关系数；h(α，τ，x)是时间尺度空间域中的滤波函数；α，τwvi(α，τ，x)是时间尺度空间域中的小波互相关原始信号。
8.如权利要求7所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，对步骤s3中时间尺度空间域的三维滤波方程进行二维傅立叶变换，得到时间频率波数域，所述时间频率波数域的三维滤波方程为：
其中，wv0(α，ω，k)是滤波结果在时间频率波数域中小波互相关经过滤波的小波互相关系数，vv(ω，k)是频率波数域中虚拟震源道集的互相关谱，p(α，ω，k)是依赖于尺度因子α的滤波函数；
9.如权利要求7所述的基于小波互相关时移地震虚拟震源波场重构方法，其特征在于，所述步骤s4中小波互相关系数的逆变换公式为：
其中，d(t)为经过处理后的地震信号，
10.一种基于小波互相关时移地震虚拟震源波场重构系统，其特征在于，包括：
上下行波分离模块，用于将垂直地震剖面vsp地震数据进行上行波和下行波分离；
小波互相关系数获取模块，用于将所述上行波和下行波都进行小波变换，选择一检波器作为虚拟震源，把所述检波器的下行波小波系数，与其它检波器的上行波小波系数互相关，得到小波互相关系数；
滤波模块，用于对所述小波互相关系数进行时间尺度域、时间尺度波数域或时间空间域滤波；
逆变换模块，用于将经过滤波的所述小波互相关系数逆变换回时空域，得到虚拟震源道集；
共中心点叠加模块，用于将所述虚拟震源道集进行共中心点叠加，得到虚拟震源激发波场。
技术总结