基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法

摘要

本发明涉及一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，首先结合矩阵的克罗内克积性质，改写3D鬼波压制算子形式，将二维大型矩阵反演问题简化为两次小型矩阵反演；然后，将预动校正条件下的三维海上拖缆资料变换至3D Radon域，估计鬼波延迟时间的搜索范围，逐道计算不同延迟时间对应峰值振幅因数，确定压制鬼波的最佳延迟时间；最后，利用构建算子进行鬼波压制。显著提高了3D Radon域鬼波压制方法的运算效率及效果；结合Radon域波场会聚性及峰值振幅因数优化策略，实现了运算效率与压制效果的平衡，同时对起伏海面的资料性有良好的适应性；最终的运算时间是常规算法的20余倍，这有利于3D鬼波压制算法在工业界的广泛应用。

权利要求

1.一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，其特征在于，包括以下步骤：a、输入地震资料、简单观测系统信息以及Radon变换参数；

b、对数据沿inline、crossline方向进行预动校正处理；

c、利用频率域三维快速抛物线Radon变换将地震资料变换至抛物线Radon域：式中，M和D代表抛物线Radon域数据的频率切片矩阵和时空域数据的频率切片矩阵，和 代表inline、crossline方向Radon变换算子， 和 为对应算子的共轭转置；

d、确定鬼波延迟时间的搜索范围，计算不同延迟时间求得的伪分离记录对应的峰值振幅因数，式中，mi(qx,qy,t)表示第i个鬼波延迟时间求得的一次波数据，rme表示mi(qx,qy,t)的均方差,qx、qy、t表示抛物线Radon域数据的坐标；

e、利用峰值振幅因数最大值对应的鬼波延迟时间Δop构建鬼波压制算子并将其变换回时空域，使用式中，Ddeghost表示时空域一次波记录的频率切片，ω表示角频率， 表示鬼波压制算子；

f、对时空域一次波记录进行反动校正，得到最终拓频结果。

说明书

基于三维快速Radon变换的海上拖缆资料自适应鬼波压制

方法

技术领域

[0001] 本发明属于地球物理勘探技术领域，具体涉及一种对采集到的海上地震资料进行频带拓宽的方法，特别涉及一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，能够对采集到的三维海上资料进行高效鬼波压制。

背景技术

[0002] 在海洋地震勘探中，有效信号的频带受到鬼波的影响会出现陷波现象，即某一频率能量剧烈降低。陷波现象导致海上地震资料频带宽度变窄，降低地震资料的纵向分辨率。鬼波压制方法是海上资料宽频处理的一大类方法，其在不同变换域计算鬼波延迟时间，根据一次波与鬼波的传播关系，将鬼波从地震资料中分离。

[0003] 目前，二维鬼波压制方法在工业界已经较为成熟，Radon域鬼波压制方法是工业界应用最成功的方法。该方法将地震资料看作不同倾角平面波子集的叠加，在不同的子集中分别根据震源、检波器沉放深度及偏移距等信息，估计一次波与鬼波的延迟时间，达到鬼波压制的目的。但该方法局限在二维条件下，实际采集的地震资料为三维空间内波场的二维切片，常规二维算法无法准确的刻画一次波与鬼波的传播关系，因此发展三维鬼波压制方法是宽频处理的重中之重。

[0004] Sun等人(2018.3D Receiver deghosting of seismic streamer data using L1 inversion and redundant extended Radon dictionary.GEOPHYSICAL PROSPECTING 66(5).987-1003)通过构建超大型稀疏矩阵将二维Radon域鬼波压制方法拓展到三维，达到三维鬼波压制的目的。然而，该算法构建的鬼波算子维度为常规算法的几万倍，大大降低了运算效率；同时，该算法要求震源、检波器沉放深度已知，但在实际采集中，沉放深度受海水等条件影响无法精确获得。因此，该算法难以满足工业生产要求。

发明内容

[0005] 本发明的目的是针对上述现有三维鬼波压制方法的不足，提供一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，对三维Radon域鬼波压制方法进行了大量的针对性改进。

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

[0007] 本发明首先对时间-空间域数据沿inline、crossline方向进行预动校正处理，引入三维快速抛物线Radon变换将地震资料变换至抛物Radon域；接着，利用粗略的观测系统信息确定鬼波延迟时间的搜索范围；将范围内的每个元素代入鬼波压制算子得到若干个伪分离记录，计算其峰值振幅因数；然后，提取峰值振幅因数最大记录对应的鬼波延迟时间构建鬼波压制算子，利用更新后的算子压制鬼波；最后，数据反变换至时间-空间域并进行反动校正，得到最终处理结果。

[0008] 一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，包括以下步骤：

[0009] a、输入地震资料、简单观测系统信息以及Radon变换参数；

[0010] b、对数据沿inline、crossline方向进行预动校正处理；

[0011] c、利用三维快速抛物线Radon变换将地震资料变换至抛物线Radon域，使用频率域三维Radon变换：

[0012]

[0013] 式中，M和D代表抛物线Radon域数据的频率切片矩阵和时空域数据的频率切片矩阵， 和 代表inline、crossline方向Radon变换算子， 和 为对应算子的共轭转置；

[0014] d、确定鬼波延迟时间的搜索范围，计算不同延迟时间求得的伪分离记录对应的峰值振幅因数，

[0015]

[0016] 式中，mi(qx,qy,t)表示第i个鬼波延迟时间求得的一次波数据，rme表示mi(qx,qy,t)的均方差,qx、qy、t表示Radon域数据的坐标；

[0017] e、利用峰值振幅因数最大值对应的鬼波延迟时间Δop构建鬼波压制算子并将其变换回时空域，使用

[0018]

[0019] 式中，表示时空域一次波记录的频率切片，表示角频率；

[0020] f、对时空域一次波记录进行反动校正，得到最终拓频结果。

[0021] 与现有技术相比，本发明的有益效果在于：

[0022] 本发明将三维快速Radon变换引入到宽频处理当中，对该方法在3D鬼波压制中的应用进行了大量的针对性改进，结合Radon域波场会聚性及峰值振幅因数优化策略，克服了常规方法在时间成本高及要求先验信息等限制，在保证鬼波压制效果的前提下，显著提高了三维Radon域鬼波压制方法的运算效率。本发明有以下特点：

[0023] 1、三维快速Radon变换的引入无需构建超大型鬼波压制矩阵，大大降低了正、反Radon变换的计算耗时，这有利于3D鬼波压制算法在工业界的广泛应用；

[0024] 2、在抛物线Radon域使用峰值振幅因数优化策略确定鬼波延迟时间，缩小了延迟时间的搜索范围，在鬼波压制的层面上实现了求解精度和求解效率的平衡；

[0025] 3、不同速度、时间的有效信息在抛物线Radon域内具有良好的分离性，三维Radon域鬼波压制方法对复杂资料具有良好的适应性。

附图说明

[0026] 图1本发明基于三维快速Radon变换的海上资料自适应鬼波压制方法流程图；

[0027] 图2a实施例1数据时空域时间切片图；

[0028] 图2b实施例1数据抛物线Radon域时间切片图；

[0029] 图3实施例1某道峰值振幅因数效果示意图；

[0030] 图4a实施例1三维鬼波压制前展示图；

[0031] 图4b实施例1三维鬼波压制前展示图；

[0032] 图5实施例1鬼波压制前后振幅谱对比图；

[0033] 图6实施例1运算效率对比图；

[0034] 图7a实施例2三维鬼波压制前展示图；

[0035] 图7b实施例2三维鬼波压制前展示图；

[0036] 图8实施例2鬼波压制前后振幅谱对比图；

[0037] 图9实施例2运算效率对比图。

具体实施方式

[0038] 下面结合附图和实例对本发明进一步的详细说明：

[0039] 本发明基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法是通过MATLAB平台实现的，已开发Windows和Linux两种环境下的程序包，地球物理数据处理人员可利用程序包进行拓频处理，具有良好的推广价值。

[0040] 一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，包括以下步骤：

[0041] a、输入地震资料、简单观测系统信息以及Radon变换参数，观测系统信息(震源、检波器沉放位置等)为确定鬼波延迟时间搜索范围的标准。

[0042] b、对数据沿inline、crossline方向进行预动校正处理。

[0043] 二维条件下，地震资料同相轴的运动学特征表现为双曲线。双曲线为二元函数，二维时空域资料仅包含inline方向信息，无法精确拟合。因此需对资料进行预动校正处理，使其运动学特征变为类抛物线，利用抛物线Radon变换进行拟合。以此类推，从稀疏表达的角度来考虑，进行inline、crossline方向的动校正处理。

[0044] 动校正速度的选取应兼顾正反校正的保幅性及稀疏表达的精度。过大时，空间方向校正量不足，同相轴在抛物线Radon域会聚性差；过小时，浅层同相轴动校正畸变严重，反动校正后存在损伤远偏移距信息的风险。经过测试，推荐使用深层同相轴均方根速度的80-90％作为全局动校正速度。

[0045] c、利用三维快速抛物线Radon变换将地震资料变换至抛物线Radon域，使用频率域三维Radon变换：

[0046]

[0047] 式中，M和D代表M和D代表抛物线Radon域数据的频率切片矩阵和时空域数据的频率切片矩阵， 和 代表inline、crossline方向Radon变换算子， 和 为对应算子的共轭转置。

[0048] 三维抛物线Radon的简化是基于矩阵克罗内克积的性质，其核心思想为将一次大型矩阵的乘法运算转换为一次小型矩阵的左乘和右乘。矩阵维度的缩小使运算效率大大提升。

[0049] d、确定鬼波延迟时间的搜索范围，计算不同延迟时间求得的伪分离记录对应的峰值振幅因数，

[0050]

[0051] 式中，mi(qx,qy,t)表示第i个鬼波延迟时间求得的一次波数据，rme表示mi(qx,qy,t)的均方差，qx、qy、t表示Radon域数据的坐标。

[0052] 相比于时空域波场，抛物线Radon域的波场具有更好的汇聚性(如图2所示)，可将其看作近似一维条件，每道记录的理论鬼波延迟时间为检波器到水面的双城走时(检波器沉放深度的二倍/海水速度)，且其最优延迟时间在空间上变化较小。因此，鬼波延迟时间的搜索范围一般为理论鬼波延迟时间的80-120％。

[0053] 计算多道伪分离记录的方式一般是以Radon域两射线参数作为主、次关键字对记录进行排序，逐道代入不同延迟时间，得到多个伪分离记录。当延迟时间准确时，鬼波能得到很好的压制；但当延迟时间不准确时，鬼波压制算子使分离记录中出现震荡。因此，利用峰值振幅因数刻画分离记录的离散程度，确定合适精确的鬼波延迟时间(例如图3所示，因数值越大表示分离记录震荡越小，最大值表示最优鬼波延迟时间)。

[0054] e、利用峰值振幅因数最大值对应的鬼波延迟时间Δop构建鬼波压制算子并将其变换回时空域，使用

[0055] 式中，Ddeghost表示时空域一次波记录的频率切片，ω表示角频率。

[0056] 确定最优鬼波延迟时间的方式是逐道进行的。首先，计算第1个延迟时间对应的伪分离记录及其峰值振幅因数，将峰值振幅因数值和对应延迟时间保留；接着，计算下一个伪分离记录，得到其峰值振幅因数，与保留的峰值振幅因数值对比，保留大值对应的延迟时间；然后，以此类推，循环延迟时间，得到单道最优延迟时间，再逐道进行循环，得到全数据最优延迟时间的二维数组；最后，对所得数组进行圆滑以提高稳定性和连续性。在圆滑时，采用以下策略：利用中值滤波剔除数组内的突变点；再利用均值滤波平滑数组。

[0057] f、对时空域一次波记录进行反动校正，得到最终拓频结果。

[0058] 实施例1

[0059] 一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，包括以下步骤：

[0060] a、输入地震资料、简单观测系统信息以及Radon变换参数；利用数值实验以对本方法的效果进行验证。震源沉放深度5m,检波器沉放深度为10m；时间采样间隔为2ms，inline方向采样间隔为12.5m,crossline方向采样间隔为10m；该数据共有inline炮集101个，crossline道集260个，时间采样点750个(图4a)。根据观测系统估算理论鬼波延迟时间Δest为67ms，并设置Radon变换参数范围qx∈[0.3,3.3],qy∈[-0.2,0.8]；

[0061] b、对数据沿inline、crossline方向进行预动校正处理，根据深层同相轴速度，选择动校正速度2250m/s；

[0062] c、利用三维快速抛物线Radon变换将地震资料变换至抛物线Radon域，使用频率域三维Radon变换将地震资料(图2a)变换至抛物线Radon域(图2b)

[0063]

[0064] 式中，M和D代表抛物线Radon域数据的频率切片矩阵和时空域数据的频率切片矩阵， 和 代表inline、crossline方向Radon变换算子， 和 为对应算子的共轭转置；

[0065] d、确定鬼波延迟时间的搜索范围Δ∈[0.9Δest,1.1Δest]，计算不同延迟时间求得的伪分离记录对应的峰值振幅因数(图3)

[0066]

[0067] 式中，mi(qx,qy,t)表示第i个鬼波延迟时间求得的一次波数据，rme表示mi(qx,qy,t)的均方差,qx、qy、t表示Radon域数据的坐标；

[0068] e、利用峰值振幅因数最大值对应的鬼波延迟时间Δop构建鬼波压制算子并将其变换回时空域，使用

[0069]

[0070] 式中，Ddeghost表示时空域一次波记录的频率切片，ω表示角频率；

[0071] f、对时空域一次波记录进行反动校正，得到最终拓频结果(图4b)。

[0072] 利用数值实验以对本方法的效果进行验证。震源沉放深度5m，检波器沉放深度为10m；时间采样间隔为2ms，inline方向采样间隔为12.5m,crossline方向采样间隔为10m；该数据共有inline炮集101个，crossline道集260个，时间采样点750个(图4a)。根据观测系统估算理论鬼波延迟时间Δest为67ms，并设置Radon变换参数范围qx∈[0.3,3.3],qy∈[-0.2,
0.8]；进行inline、crossline方向预动校正处理，根据深层同相轴速度，选择动校正速度
2250m/s；将地震资料变换至抛物线Radon域，逐道更新最优鬼波延迟时间，构建鬼波压制算子，压制Radon域鬼波能量，对比压制前后的振幅谱与计算时间，评估本发明所提方法的有效性。完整的流程与简要的结果分析如下：

[0073] 设置初始参数，进行预动校正处理，使用三维快速抛物线Radon变换将地震资料(图2a)变换至抛物线Radon域(图2b)；

[0074] 确定鬼波延迟时间的搜索范围Δ∈[0.9Δest,1.1Δest]，计算不同延迟时间对应的峰值振幅因数(图3)，利用最优鬼波延迟时间构建鬼波压制算子，将数据变回时空域，进行反动校正得到最终鬼波压制结果(图4)，同相轴子波旁瓣能量减弱，鬼波得到压制，但炮集中不够明显；

[0075] 计算inline方向前50个道集鬼波压制前后的平均振幅谱(图5，图中，虚线表示鬼波压制前，实线表示鬼波压制后)，鬼波导致的“锯齿状”陷波点得到良好的补偿，压制后频带明显拓宽，证明了方法的有效性；记录常规方法与快速自适应鬼波压制方法的运算时间(图6，其中左侧的分量1为常规三维Radon域鬼波压制方法计算1炮的计算时间；右侧的分量2为本发明提出的三维Radon域鬼波压制方法计算1炮的计算时间)，快速自适应鬼波压制方法优势明显，提升运算效率近30倍。

[0076] 实施例2

[0077] 一种基于三维快速Radon变换的海上拖缆资料自适应鬼波压制方法，包括以下步骤：

[0078] a、输入地震资料、简单观测系统信息以及Radon变换参数；利用实际资料实验以对本方法的效果进行验证；该资料震源沉放深度10m,检波器沉放深度为20m；时间采样间隔为4ms，inline方向采样间隔为12.5m,crossline方向采样间隔为20m；该数据共有inline炮集
48个，crossline道集120个，时间采样点500个(图7a)。根据观测系统估算理论鬼波延迟时间Δest为135ms，并设置Radon变换参数范围qx∈[0.3,3.3],qy∈[-0.1,0.4]；

[0079] b、对数据沿inline、crossline方向进行预动校正处理，根据深层同相轴速度，选择动校正速度2500m/s；

[0080] c、利用三维快速抛物线Radon变换将地震资料变换至抛物线Radon域，使用频率域三维Radon变换将地震资料变换至抛物线Radon域

[0081]

[0082] 式中，M和D代表M和D代表抛物线Radon域数据的频率切片矩阵和时空域数据的频率切片矩阵， 和 代表inline、crossline方向Radon变换算子， 和 为对应算子的共轭转置；

[0083] d、确定鬼波延迟时间的搜索范围Δ∈[0.8Δest,1.2Δest],计算不同延迟时间求得的伪分离记录对应的峰值振幅因数

[0084]

[0085] 式中，mi(qx,qy,t)表示第i个鬼波延迟时间求得的一次波数据，rme表示mi(qx,qy,t)的均方差,qx、qy、t表示Radon域数据的坐标；

[0086] e、利用峰值振幅因数最大值对应的鬼波延迟时间Δop构建鬼波压制算子并将其变换回时空域，使用

[0087]

[0088] 式中，Ddeghost表示时空域一次波记录的频率切片，ω表示角频率；

[0089] f、对时空域一次波记录进行反动校正，得到最终拓频结果(图7b)。

[0090] 对实际地震资料进行快速自适应压制鬼波处理，该资料震源沉放深度10m，检波器沉放深度为20m；时间采样间隔为4ms，inline方向采样间隔为12.5m,crossline方向采样间隔为20m；该数据共有inline炮集48个，crossline道集120个，时间采样点500个(图7a)。采用与上一算例相近的参数：理论鬼波延迟时间Δest为135ms，Radon变换参数范围qx∈[0.3,3.3],qy∈[-0.1,0.4]，预动校正速度2500m/s，测试本发明中的方法对实际数据的有效性。

[0091] 在频带拓宽方面，提取三维数据体的不同切片(图7)，同相轴旁瓣能量减弱(crossline:-500m,inline:0m,time:1.8s)，验证了算法的有效性；压制前后平均振幅谱的对比图(图8，图中，虚线表示鬼波压制前，实线表示鬼波压制后)同样证明了算法的有效性，然后，压制后的振幅谱(红色实线)在20-50Hz的位置仍存在一定的陷波现象；产生算子的原因是实际资料在传播过程因地层吸收，子波主频随传播时间增加而减少，鬼波压制算子在高频分量中不适用。在计算效率方面，由图9可看出，快速自适应方法提升运算效率近20倍(图中，左侧分量1为常规三维Radon域鬼波压制方法计算1炮的计算时间；右侧分量2为本发明提出的三维Radon域鬼波压制方法计算1炮的计算时间)。

[0092] 结合以上两个点，在本算例条件下，三维快速自适应鬼波压制技术在实际地震资料处理中，兼顾了拓宽频带和降低运算成本两方面，推广意义大。