一种特征波域的浅水多次波模型指向性预测方法

摘要

本发明涉及一种特征波域的浅水多次波模型指向性预测方法，属于地球物理勘探技术领域，解决了现有浅水多次波预测方法预测结果存在假象干扰的问题。包括以下步骤：在特征波域对地震数据进行分解，得到特征波数据；提取特征波传播方向预测算子；根据上述提取的特征波传播方向预测算子进行特征波域指向性预测。该方法通过引入特征波分解与合成，利用了地震数据的方向信息，使得多次波预测时能够避免虚假路径，减少了假象噪音干扰，提高了预测精度；同时，预测的浅水多次波形态、幅值都与实际相符，有效提升了目标靶区地震资料的信噪比和分辨率，另外，本方法可应用于三维地震资料浅水多次波的压制,大大提升了浅水区油气勘探的数据处理效果。

说明书

技术领域

本发明涉及地球物理勘探技术领域，尤其涉及一种特征波域的浅水多次波模型指向性预测方法。

背景技术

在近海地震勘探中，浅水区多次波特别发育，严重影响地震资料的可解释性，目前，压制多次波的常规方法主要有：基于一次波和多次波空间差异特性的滤波类方法(如：线性/非线性Radon变换压制多次波等)，波动理论多次波压制方法(如：预测反褶积、数据驱动的波动理论预测-减去类方法、模型驱动的波动理论预测-减去类方法等)。

由于浅水近偏移距数据缺失、信噪比低、水体相关多次波异常发育、多次波与有效波严重重合等原因，传统压制多次波的方法对于浅水多次波都存在着缺点和不足：如由于数据的信噪比与完备性，数据驱动类方法(如SRME)不能有效地预测水体相关多次波；模型驱动的波动理论预测-减去类方法：在水体模型作为先验信息的基础上，根据波动理论来预测水体相关多次波。此类方法通过波动方程外推来实现多次波的预测，弱化了数据驱动类方法对数据完整性的要求，对浅水多次波压制有一定的适应性。

但是，现有的模型驱动类多次波压制方法由于没有考虑方向信息，虽然对多次波贡献道集叠加孔径进行了优化，预测结果依旧含有大量的假象噪音干扰；同时，现有的方法一般仅能预测多次波的运动学信息，预测的多次波的振幅和相位等动力学信息与真实多次波存在差异，为之后的预测减去带来干扰，影响了多次波的压制效果。

发明内容

鉴于上述的分析，本发明旨在提供一种特征波域的浅水多次波模型指向性预测方法，用以解决现有浅水多次波预测方法预测结果存在假象干扰的问题。

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

提供了一种特征波域的浅水多次波模型指向性预测方法，包括以下步骤：

在特征波域对地震数据进行分解，得到特征波数据；

提取特征波传播方向预测算子；

根据上述提取的特征波传播方向预测算子进行特征波域指向性预测。

本发明有益效果如下：

本发明方法针对浅水区域的特点，将多次波预测问题转换为两个核心问题，一个是特征波变换，一个是多次波预测问题，通过引入特征波分解与合成，利用了地震数据的方向信息，使得多次波预测时能够避免虚假路径，减少了假象噪音干扰，提高了预测精度。同时，预测的浅水多次波形态、幅值都与实际相符，为下一步的多次波匹配减去奠定了基础，进而有效提升了目标靶区地震资料的信噪比和分辨率，为高分辨率地震成像提供了支持。另外，本方法可应用于三维地震资料浅水多次波的压制，大大提升了浅水区油气勘探的数据处理效果，具有重要的实践意义。

在上述方案的基础上，本发明还做了如下改进：

进一步，所述根据上述提取的特征波传播方向预测算子进行特征波域指向性预测，包括：

根据每一特征波传播方向预测算子生成每个局部空间窗的特征波域多次波预测模型；

将上述所有局部空间窗的特征波域多次波预测模型叠加，得到全区域特征波域多次波预测模型；

将上述全区域特征波域多次波预测模型变换回到时空域，得到最终的多次波预测模型。

进一步，所述提取特征波传播方向预测算子，包括：

将时空域的地震数据划分为若干个局部空间窗；

选取每个局部空间窗所对应的特征波数据；

按上述每个特征波数据的方向在水体模型内做单点射线追踪，得到该局部空间窗的特征波传播方向预测算子。

进一步，还包括：对局部空间窗的数据通过加权系数进行约束，所述加权系数选用高斯函数窗，在中心射线路径上数据权重＝1，偏离中心射线的数据权重递减。

进一步，所述局部空间窗的特征波域多次波预测模型，公式为：

式中，x_r为待预测多次波所在检波点的位置，x′为二次源的位置，p为射线参数，j为射线参数的个数，x₀为参考道的位置，s’为二次源所在平面，ω为频率，U为平面波源，G为格林函数，n为外法线方向。

进一步，对于每个局部空间窗所对应的单频ω，单方向p_ip的特征波数据，存在以下预测关系：

式中，m(x_s，x_n，p_ip，ω)表示为炮点在x_s时，检波点在x_i时，射线参数为p_ip的特征波域的多次波预测模型；s(x_s，x′_n，p_ip，ω)表示为炮点在x_s时，位于x′_i的各检波点接收到的数据；g(x′_i，x_j，p_ip，ω)为预测算子矩阵，即从二次源x′_i传至检波点x_j的格林函数。

进一步，所述在特征波域对地震数据进行分解，包括选取特征波基向量组：

A_(k，l)＝exp(iωp_l(x_k-x₀))，

其中，p_l为射线参数，l为射线参数的个数，x_k为地震数据的位置，k为地震数据的道数，x₀为参考道的位置。

进一步，所述在特征波域对地震数据进行分解，还包括：采用预条件共轭梯度法求解特征波域系数，公式为：

式中，d是地震数据向量，s是待求的特征波域系数，A是所选取的特征波基向量组；是加权矩阵的近似，为的逆矩阵，m为中间参数。

进一步，还包括：对原始地震数据进行去噪处理及数据重构。

进一步，所述地震数据为浅水区数据，所述多次波包括：鬼波、水体相关多次波、自由表面相关多次波。

本发明中，上述各技术方案之间还可以相互组合，以实现更多的优选组合方案。本发明的其他特征和优点将在随后的说明书中阐述，并且，部分优点可从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过说明书、权利要求书以及附图中所特别指出的内容中来实现和获得。

附图说明

附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。

图1为本方法实施例中特征波域的浅水多次波模型指向性预测方法流程图；

图2为常规模型驱动类多次波预测方法示意图；

图3为本发明实施例中特征波指向性预测方法示意图；

图4为波场传递路径衔接示意图；

图5为实际OBC观测系统中的水体多次波传播路径示意图；

图6为理论OBC观测系统的单炮数据；

图7为采用常规方法对图6中的数据进行多次波预测的结果示意图；

图8为采用实施例1中的方法对图6中的数据进行多次波预测结果示意图；

图9为实际OBC观测系统中的单炮数据；

图10为采用实施例1中方法对图9中的数据进行多次波预测的结果示意图；

图11为对图9中数据中的多次波压制后的结果示意图。

具体实施方式

下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。

实施例1

本发明的一个具体实施例，公开了一种特征波域的浅水多次波模型指向性预测方法，如图1所示，包括以下步骤：

步骤S1、在特征波域对地震数据进行分解，得到特征波数据；

步骤S2、提取特征波传播方向预测算子；

步骤S3、根据上述提取的特征波传播预测算子进行特征波域指向性预测。

与现有技术相比，本实施例提供的特征波域的浅水多次波模型指向性预测方法，针对浅水区域的特点，将多次波预测问题转换为两个核心问题，一个是特征波变换，一个是多次波预测问题，通过引入特征波分解与合成，利用了地震数据的方向信息，使得多次波预测时能够避免虚假路径，减少了假象噪音干扰，提高了预测精度。同时，预测的浅水多次波形态、幅值都与实际相符，为下一步的多次波匹配减去奠定了基础，进而有效提升了目标靶区地震资料的信噪比和分辨率，为高分辨率地震成像提供了支持。另外，本方法可应用于三维地震资料浅水多次波的压制，大大提升了浅水区油气勘探的数据处理效果，具有重要的实践意义。

需要说明的是，本实施例方法通过引入特征波概念(特征波：时空局部、单震相、带方向的带限波场)，成功简化了多次波与一次波的预测关系。在特征波域，多次波可被视为有效波作为二次源的波场外推结果，外推算子(传播方向预测算子)是地质模型的格林函数；数据可稀疏表达为少量几个局部平面波；预测关系由经典的二次点源球面波叠加精简为指向性更强的二次面源局部平面波叠加；浅水区的格林函数简化为射线理论下的射线走时。同时，再结合压缩感知稀疏反演技术，即使在近偏移距数据缺失时，仍能显著提高浅水区域的多次波模型预测精度。

具体来说，

在步骤S1中、在特征波域对原始地震数据进行分解，得到特征波数据；其中，原始地震数据为浅水区地震数据。

根据射线走向对地震数据进行稀疏表达，可以在特征波场域对地震数据做分解；其中，特征波场为局部时间空间域的单震相的方向带限波场。

需要说明的是，在对原始地震数据进行稀疏表达之前，可以对原始地震数据进行预处理。通过采用现有的多种去噪方法去除地震数据的噪声干扰，以便提高地震数据的分辨率，达到更好的多次波预测效果。

另外，由于浅水区地震数据存在近偏移距数据缺失及坏道的问题，可以采用基于波动方程的地震数据构建方法、基于预测滤波的地震数据构建方法、基于压缩感知框架的数学变换地震数据构建方法对地震数据进行重构，以便得到完整、准确的地震数据。

将时空域的地震数据转换到特征波域，通过选取一组正交基，将地震数据投影到基函数空间中，形成有限个数的、幂指数衰减的投影系数值，使得变换域中的地震数据具有稀疏性质。具体地。将转换过程可表述为如下优化问题：

Min||s||_o

其中，d是地震数据向量，s是待求的特征波域系数，ε为约束系数，A是所选取的特征波基向量组，如下式所示：

A_(k，l)＝exP(iωp_l(x_k-x₀))，

其中，p_l(l＝0,1，…)为射线参数，l为射线参数的个数，x_k(k＝0,1…)为地震数据(或局部空间窗的数据)的位置，k为地震数据(或局部空间窗的数据)的道数，x₀为参考道的位置，优选的，选择地震数据(或局部空间窗)中心道的位置作为参考道。

求解得到的各个方向的特征波分量s(p_l)，这实际上赋予了数据额外的方向信息，相当于将数据(时空、频率)升维到更高维空间(空间、频率、方向)，即d(x，ω)→d(x，p，ω)。

为了便于求解计算，本实施例采用加权L₂范数逼近L₀范数，则将上述待求优化问题近似于如下无约束优化问题：

其中，d是地震数据向量，s是待求的特征波域系数，A是所选取的特征波基向量组；P是加权矩阵，可以取为一个对角阵，对角元是待求的特征波域系数s中各分量能量的倒数，如下式所示：

其中，s_n为待求的特征波域系数s中第n个分量。

由于s本身为待求系数，本实施例通过近似表示得到，即：

其中，为待求的特征波域系数s中第n个分量的近似值。

则上述待求无约束优化问题转化为：

式中，d是地震数据向量，s是待求的特征波域系数，A是所选取的特征波基向量组；是加权矩阵的近似。

进一步，可采用共轭梯度等算法直接对上述转化后的待求无约束优化问题求解其法方程，得到待求的特征波域系数s。

为了提高收敛速度，本实施例通过变量代换，采用预条件共轭梯度法进行求解：

先用预条件共轭梯度法求解出中间参数m，再由导出待求的特征波域系数s。需要说明的是，求取的s即为地震数据的稀疏表达，其的物理意义为地震数据中各个方向的平面波的复振幅。

步骤S2、提取特征波传播方向预测算子。对每个特征波分量，按其方向做射线追踪，得到精确的符合物理的指向性预测算子。

具体地，将时空域的地震数据划分为若干个局部空间窗；取出第k个局部空间窗所对应的特征波数据，按其方向在水体模型内做单点射线追踪使其在水体中发生反射并再次传至检波面，得到射线意义下的特征波传播预测算子。在特征波域，数据有额外的方向信息，故能按照其真实方向外推，并在折射、反射界面处严格遵守snell定律，从而避免了虚假波路径的产生，提高了预测质量。

需要说明的是，局部空间窗的尺寸可以根据实际数据的大小、精度要求、硬件处理能力等进行设定，优选的，选取时间维度为500ms，空间维度300米半径的窗口。

为了使预测的多次波的动力学信息(相位、振幅)更加符合真实情况，本实施例通过对局部空间窗的数据通过加权系数进行约束，加权系数用于表征每个特征波域二次源对多次波贡献道集的贡献大小。在菲涅尔带内的真实二次源处，加权系数较大；而菲涅尔带之外的其它非真实二次源处，权系数迅速衰减。优选的，加权系数选用高斯函数窗，在中心射线路径上权重＝1，偏离中心射线的二次源权重递减。

进一步，针对自由表面相关多次波，还可以采用二次源处的自由表面加权算子(包含自由表面反射系数、子波的逆、检波器矩阵的共同作用的自由表面算子及上述加权系数的两者结合)对局部空间窗的数据进行约束。

在特征波域编码预测多次波的基础上，引入了特征波域加权系数及自由表面加权算子，实现了特征波域多次波运动学和动力学信息的统一预测。

步骤S3、根据上述提取的特征波传播预测算子进行特征波域指向性预测；将步骤S2中提取的第k个局部空间窗的传播预测算子作用于对应的特征波域数据，进行特征波域指向性预测，得到该局部空间窗所产生的多次波模型，进而得到全部区域的多次波预测模型。

具体地，步骤S301、根据特征波传播方向预测算子生成每个局部空间窗的特征波域多次波预测模型；

将局部空间窗的特征波传播方向预测算子作用于对应的特征波域的数据，生成二次源多次波贡献道集，通过所有二次源波场多次波贡献道集的叠加，得到该局部空间窗的特征波域多次波预测模型。

现有模型驱动压制多次波方法，如图2所示，将平面上每一点都作为二次震源(点状偶极源)，激发球面波传播到平面后的全空间。在平面上每个可能的向下反射点处将输入记录与预测算子褶积，形成多次波贡献道集，通过多次波贡献道集的叠加来预测多次波。只有对应于真实向下反射点的多次波贡献道集叠加才会产生相长干涉，其它部分的多次波贡献道集叠加则会产生相消干涉，影响多次波预测的精度，而真实的向下反射点位置或多次波贡献道集叠加孔径均是未知的，虽然对叠加孔径进行了优化，但是叠加后的数据依旧存在大量假象噪音干扰。

本实施例中，如图3所示，平面上每一点都作为二次源(有方向的局部平面波源)传播到平面后的射线路径及其邻域；以特征波域的观测波场和满足斯奈尔定律的一次反射格林函数作为二次源波场和特征波传播方向预测算子，进而通过选取符合斯奈尔定律的二次源贡献形成多次波贡献道集，实现在特征波域对多次波贡献道集的自动优化，以去除假象干扰，提高多次波的预测精度。

对于特征波域中单频ω，单方向p_ip的数据，有如下预测关系

式中，(m(x_s，x_n，p_ip，ω))表示为炮点在x_s时，检波点在x_i(i＝0,1,2…n)时，射线参数为p_ip的特征波域的多次波预测结果。(s(x_s，x′_n，p_ip，ω))表示为炮点在x_s时，位于x′_i(i＝0,1,2…n)的各检波点接收到的数据。G(g(x′_i，x_j，-p_ip，ω))为预测算子矩阵，即从二次源x′_i传至检波点x_j的格林函数。

上式描述了符合snell定律的波路径连接过程，如图4所示，保证了在反射面元(二次源)x′_i处，入射角等于反射角，有效避免了虚假反射。

具体地，波动理论Rayleigh积分可表示为：

式中，x_r为待预测多次波所在道的位置，s’为积分球面，x′为二次源的位置。

在特征波域对二次源进行展开：

式中，x′为二次源的位置，p为射线参数，j为射线参数的个数，x₀为参考道的位置。

将上述“特征波展开”带入Rayleigh积分：

式中，x_r为待预测多次波所在道的位置，x₀为参考道的位置，x′为二次源的位置，p为射线参数，s’为二次源所在平面。

交换积分、求和顺序，得到特征波域多次波预测公式：

式中，x_r为待预测多次波所在检波点的位置，x′为二次源的位置，p为射线参数，j为射线参数的个数，x₀为参考道的位置，s’为二次源所在平面，ω为频率，U为平面波源，G为格林函数，n为外法线方向。

有上式可知，可先计算单个平面波成分的Rayleigh积分，再将各个平面波成分累加，得到特征波域多次波预测模型。

对于单方向、单个频率的平面波，其Rayleigh积分的结果有解析解—平面波外推(均匀介质非反射情形)：

式中，x_r为待预测多次波所在道的位置，x₀为参考道的位置，x′为二次源的位置，p0为射线参数。

需要说明的是，在求取特征波域多次波预测模型过程中，特征波域的多次波贡献道集可以表示为：

式中，x_r为待预测多次波所在道的位置，x′为二次源的位置，p为射线参数，j为射线参数的个数。需要说明的是，得到的多次波贡献道集相当于忽略了不符合反射定律的所有射线路径之贡献。

步骤S302、输出多次波模型预测结果，将上述所有局部空间窗的特征波域多次波预测模型叠加，得到全区域特征波域多次波预测模型；

具体地，对地震数据中所有的局部空间窗进行指向性预测处理，并将各个局部空间窗预测得到的多次波模型叠加得到最终的浅水多次波预测结果。

需要说明的是，可以将各个局部空间窗的特征波域的多次波模型叠加，得到全区域特征波域多次波预测模型，并将该叠加后特征波域多次波模型变换回到时空域，得到时空域预测的多次波模型，即多次波模型最终预测结果；还可以，将各个局部空间窗的特征波域的多次波模型分别变换到时空域，在时空域进行叠加，得到最终的多次波预测模型。

实施例2

实施例1中指向性预测方法在海上地震勘探中有重要价值；具体来说，在浅水区域的OBC、OBN、拖缆数据的高分辨率处理中，本方法有助于精确压制浅水多次波，提高数据质量。其中，压制的浅水多次波主要主要包括鬼波、水体相关多次波和自由表面相关多次波；可以对上述浅水区的多次波进行单独压制，也可以将自由表面相关多次波和鬼波统一进行压制。

本实施例将实施例1中的浅水多次波模型指向性预测方法应用到OBC观测数据中，用于说明该方法的实际应用价值及有益效果。在OBC观测系统中，检端水体相关多次波等价于检鬼波，源端水体相关多次波与常规拖缆观测系统中的源端水体相关多次波相同。此时，二次源波场是位于海底的观测波场，射线路径如图5中实线所示；预测方向传播算子为描述波场由二次源位置向上传播至自由表面再向下反射至位于海底的检波点处的格林函数，射线路径如图5中虚线所示。

本实施例中通过实施例1中的方法，对OBC观测系统理论数据及实际数据进行处理，并与常规浅水多次波压制方法进行对比，如图6-8所示。能明显看出本方法预测得到的多次波含有更少假象，有利于提高地震数据的质量。

另外，采用实施例1中的方法，对三维浅水OBC实际观测数据进行处理，特征波域的指向性预测方法可以较准确地预测出水体相关多次波，证明了本发明方法的有效性，如图9-10所示。接下来，将预测结果从原始记录中自适应减去，得到了水体相关多次波的压制结果，如图11所示。从图11中的箭头所示区域对比说明，本发明方法对浅水(水体相关)多次波的压制有良好效果。

本领域技术人员可以理解，实现上述实施例方法的全部或部分流程，可以通过计算机程序指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。