CN114609671A - 鬼波衰减方法、装置、计算机设备及可读存储介质 - Google Patents

鬼波衰减方法、装置、计算机设备及可读存储介质 Download PDF

Info

Publication number
CN114609671A
CN114609671A CN202011420878.2A CN202011420878A CN114609671A CN 114609671 A CN114609671 A CN 114609671A CN 202011420878 A CN202011420878 A CN 202011420878A CN 114609671 A CN114609671 A CN 114609671A
Authority
CN
China
Prior art keywords
matrix
ghost
regularization parameter
iteration
wave
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202011420878.2A
Other languages
English (en)
Inventor
王兆旗
范国章
李立胜
李林
王彬
杨涛涛
叶月明
杨志力
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN202011420878.2A priority Critical patent/CN114609671A/zh
Publication of CN114609671A publication Critical patent/CN114609671A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Oceanography (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明实施例提供了一种鬼波衰减方法、装置、计算机设备及可读存储介质,其中,该方法包括:获取海上水平拖缆地震资料;基于平面波传播原理,根据海上水平拖缆地震资料建立τ‑p域拖缆观测的总波场与海水表面上行波波场之间的关系方程;利用最小平方残差方法反演求解关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。该方案实现了采用混合最小平方残差方法来反演求解关系方程,有利于使得迭代解逼近真实解,有利于提高反演求解结果的准确性,有利于提高计算效率的同时改善鬼波的压制效果。

Description

鬼波衰减方法、装置、计算机设备及可读存储介质
技术领域
本发明涉及海洋地震勘探技术领域,特别涉及一种鬼波衰减方法、装置、计算机设备及可读存储介质。
背景技术
高保真、宽频带和准确成像的地震资料处理技术是高品质地震成像的基础,针对性处理措施是提高地震资料处理质量的核心。在在海洋地震勘探中,鬼波的压制是一个关键,鬼波不仅干扰有效波的成像,而且减宽了原始地震信号的频带,给地质体的有效识别造成了严重困难。近几年,鬼波压制问题日益受到业界重视,已成为海上地震资料宽频处理中的核心技术之一。
目前,国内外的地球物理学家分别从海上地震采集方式和地震数据处理方法两方面进行了深入的研究,一是在野外采集中采用合适的震源组合、电缆形态或电缆组合减小鬼波的影响;二是在资料处理过程中采用合适的算法从实测数据中消除鬼波。众多学者开展了鬼波压制处理算法研究,现有技术中提出了利用常规偏移和镜像偏移联合反褶积的方法衰减鬼波,现有技术中还通过反演得到鬼波最优延迟时间,在f-x-y域和τ-p域实现了鬼波压制。现有技术中还有利用斜缆数据中鬼波与一次反射波的视速度的差异将鬼波与一次反射波分离的方法。现有技术中还有方案假设海水反射系数是频率和慢度的函数,估计鬼波的最优延迟时间进行鬼波压制。现有技术中还有方案利用最小平方残差方法(LSMR)反演求解方程得到海水表面观测的上行波波场,但是没有考虑到LSMR算法的非收敛性,仍然存在计算效率低、处理结果精度低的问题。
发明内容
本发明实施例提供了一种鬼波衰减方法,以解决现有技术中鬼波衰减存在的计算效率低、处理结果精度低的技术问题。该方法包括:
获取海上水平拖缆地震资料;
基于平面波传播原理,根据所述海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
本发明实施例还提供了一种鬼波衰减装置,以解决现有技术中鬼波衰减存在的计算效率低、处理结果精度低的技术问题。该装置包括:
海上地震资料获取模块,用于获取海上水平拖缆地震资料;
方程建立模块,用于基于平面波传播原理,根据所述海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
方程求解模块,用于利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
数据处理模块,用于对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
本发明实施例还提供了一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述任意的鬼波衰减方法,以解决现有技术中鬼波衰减存在的计算效率低、处理结果精度低的技术问题。
本发明实施例还提供了一种计算机可读存储介质,所述计算机可读存储介质存储有执行上述任意的鬼波衰减方法的计算机程序,以解决现有技术中鬼波衰减存在的计算效率低、处理结果精度低的技术问题。
在本发明实施例中,基于平面波传播原理,根据海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程后,利用最小平方残差方法反演求解关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场,进而对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。提出了在最小平方残差方法迭代过程加入正则化参数,将正则化参数应用于最小平方残差方法的子问题,实现了采用混合最小平方残差方法来反演求解关系方程,有利于使得迭代解逼近真实解,进而有利于提高反演求解结果的准确性,有利于提高计算效率的同时改善鬼波的压制效果。
附图说明
此处所说明的附图用来提供对本发明的进一步理解,构成本申请的一部分,并不构成对本发明的限定。在附图中:
图1是本发明实施例提供的一种鬼波衰减方法的流程图;
图2是本发明实施例提供的一种拖缆地震数据记录的总波场与海底地层的上行反射波场之间的关系示意图;
图3是本发明实施例提供的一种层状模型及模拟单炮记录的示意图;
图4是本发明实施例提供的一种分别采用LSMR算法和混合LSMR算法对图3中(b)单炮记录进行鬼波压制测试的结果对比图;
图5是本发明实施例提供的一种图3中(b)单炮记录频谱图及其分别采用LSMR算法和采用混合LSMR算法压制鬼波后的记录频谱图;
图6是本发明实施例提供的一种测试不同正则化参数的鬼波压制结果示意图;
图7是本发明实施例提供的一种包含鬼波的实际地震记录及鬼波压制后的实际地震记录的示意图;
图8是本发明实施例提供的一种计算机设备的结构框图;
图9是本发明实施例提供的一种鬼波衰减装置的结构框图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚明白,下面结合实施方式和附图,对本发明做进一步详细说明。在此,本发明的示意性实施方式及其说明用于解释本发明,但并不作为对本发明的限定。
在本发明实施例中,提供了一种鬼波衰减方法,如图1所示,该方法包括:
步骤102:获取海上水平拖缆地震资料;
步骤104:基于平面波传播原理,根据所述海上水平拖缆地震资料建立拉东(τ-p)域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
步骤106:利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
步骤108:对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
由图1所示的流程可知,在本发明实施例中,在本发明实施例中,基于平面波传播原理,根据海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程后,利用最小平方残差方法反演求解关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场,进而对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。提出了在最小平方残差方法迭代过程加入正则化参数,将正则化参数应用于最小平方残差方法的子问题,实现了采用混合最小平方残差方法来反演求解关系方程,有利于使得迭代解逼近真实解,进而有利于提高反演求解结果的准确性,有利于提高计算效率的同时改善鬼波的压制效果。
具体实施时,获取海上水平拖缆地震资料后,在建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程之前,可以先对海上水平拖缆地震资料开展叠前预处理,叠前预处理可以主要包括以下步骤:野外地震数据解编为内部格式,观测系统加载;坏炮与坏道编辑,枪深、缆深校正;气泡压制;初始球面扩散补偿;低频滤波;涌浪噪音衰减;线性噪音衰减。
具体实施时,基于平面波传播原理,根据所述海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程的过程可以参考现有技术实现,例如,排除船舶反射等环境背景噪声的影响因素,海洋拖缆地震数据所记录的波场a(xi,t)可以视为来自海底地层的上行反射波场u(xi,t)和该上行反射波场传播到海面后再由海面反射向下传播的下行波场d(xi,t)(即鬼波)之和:
a(xi,t)=u(xi,t)+d(xi,t) (1)
如图2所示,D点在某时刻接收到的总波场p是反射波传播到D点的上行波场u与由传播到A点再由A点反射到D点的下行波场d之和。上行波场u到达C点与上行波场u到达B点的时差
Figure BDA0002822300630000041
为:
Figure BDA0002822300630000042
上式中,θm为第m个出射角,v为海水速度,pm为第m个射线参数,zi为第i道相对于海面的沉放深度。
进一步可以得到:D点在t时刻接收到的总波场p可以视作C点在
Figure BDA0002822300630000051
时刻接收到的上行波场
Figure BDA0002822300630000052
与A点在
Figure BDA0002822300630000053
时刻接收到的上行波场
Figure BDA00028223006300000513
乘以海面反射系数之后求和,即下式:
Figure BDA0002822300630000054
上式中:xi为该检波器所在道的炮间距,R代表海水反射系数。
假设海水弹性参数在横向上均匀不变,则A、C两点所接收到的上行波u0(xi-Δxi,m,t)、u0(xi+Δxi,m,t)与B点接收到的上行波
Figure BDA0002822300630000055
分别存在如下关系:
Figure BDA0002822300630000056
Figure BDA0002822300630000057
联合式(3)和(4),结合频率域线性τ-p反变换,我们可以得到:
Figure BDA0002822300630000058
上式中:A(xit)、V(pm,ω)分别代表总波场、上行波场u0(xi,t)在频率域拉东变换的结果,ω代表频率,m代表慢度道的总个数,j=1,2,…m,
Figure BDA0002822300630000059
如图2所示。
可将上式(5)写成矩阵形式,即关系方程:A=GU
其中:
Figure BDA00028223006300000510
A、G、U分别代表总波场、拉东反变换算子以及上行波场的矩阵表示。
G为拉东反变换算子,可分解为上行波τ-p变换算子Gu和下行波τ-p变换算子Gd之和。即:
Figure BDA00028223006300000511
式(7)中,
Figure BDA00028223006300000512
代表拉东反变换算子G的分解项,pj代表第j个射线参数。
具体实施时,得到τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程之后,即可求解关系方程得到压制鬼波后的一次波,现有技术一般利用最小平方残差方法(LSMR)反演求解关系方程(6),即求解线性系统Ax=b和最小平方问题min||Ax-b||2,得到海水表面观测的上行波波场。
LSMR算法的具体实现就是通过Golub-Kahan过程实现双对角化。给定矩阵A和向量b,GK过程是将矩阵(b A)变换为上双对角形式(β1e1Bk)的迭代过程。选择两个正常数α11,使得数据空间正交基矢量uk和解空间正交基矢量vk的范数等于1。进行初始化β1=||b||2,u1=b/β1,α1v1=ATu1,GK过程的第k次迭代计算正交向量uk+1和vk+1,使得:
βk+1uk+1=Avkkuk αk+1vk+1=ATuk+1k+1vk (8)
在k步之后,我们有矩阵
Vk=(v1,v2,…,vk);Uk=(u1,u2,…,uk),并得到一个新的双对角矩阵Bk
Figure BDA0002822300630000061
对于LSMR,在每次迭代时选择yk以最小化||ATrk||2,此处rk为迭代残差,rk=b-Axk,将解决方案定义为xk=Vkyk
Figure BDA0002822300630000062
上式中:e代表单位矩阵,ek+1表示单位矩阵中第k+1列向量,y代表子问题的解决方案,yk代表第k步之后y的值。
Figure BDA0002822300630000063
αk、βk、βk+1、αk+1分别代表第k、k+1步之后α和β的值。
但是,本申请发明人发现,由于LSMR表现出半收敛行为,重建误差最初减小但在某些点开始增加,近似解在迭代初始阶段逼近真实解,但是在某一步之后就偏离真实解.并且
Figure BDA0002822300630000064
在更多迭代中变得病态,因此,提出了我们需要使用Tikhonov正则化来解决LSMR子问题,即在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,进而得到海水表面的上行波波场,例如:
Figure BDA0002822300630000071
具体实施时,提供了多种在每次迭代中确定正则化参数的方法,例如,
第一种:基于广义交叉验证(GCV)确定正则化参数:
通过期望极小化均方差
Figure BDA0002822300630000072
获取正则化参数λ;
Figure BDA0002822300630000073
上式(12)中:A=U∑VT代表A的奇异值,其中Σ=diag(σ1,σ2,...,σn)是包含A的奇异值的对角矩阵,σ1≥σ2≥···≥σn≥0,并且U的列ui和V的列vi分别包含左、右奇异向量,GA,b(λ)表示正则化参数矩阵,
Figure BDA0002822300630000074
代表矩阵A的广义逆,xλ代表模拟数据,
Figure BDA0002822300630000075
表示矩阵A和模拟数据组成的矩阵,b代表实测数据,A表示总波场的矩阵,I表示单位矩阵,n表示迭代总次数。
第二种:基于偏差原理(DP)确定正则化参数:
要求噪声||e||信息,其核心思想就是残差
Figure BDA0002822300630000076
不可能小于右端项的误差。即:
Figure BDA0002822300630000077
其中,δ是用户定义的安全因子,其值略大于1。
第三种:基于无偏预测风险估计(UPRE)确定正则化参数:
其思想是选择使预测风险的期望值的无偏估计最小化的正则化参数,即
对于UPRE,其思想是选择正则化参数来最小化预测风险的期望值的无偏估计,即通过以下公式(14)确定正则化参数:
Figure BDA0002822300630000078
上式中,σ2表示噪声方差,UA,b(λ)表示待求的正则化参数矩阵。
具体实施时,为了进一步提高计算效率,在本实施例中,还提出了在最小平方残差方法迭代过程中,确定终止迭代次数,当迭代次数达到终止迭代次数时,终止迭代。
具体的,通过最小化
Figure BDA0002822300630000079
得到以下公式(即GCV函数)来确定终止迭代次数k:
Figure BDA0002822300630000081
上式中:λ表示正则化参数,I表示单位矩阵,n表示迭代总次数,G(k)代表终止迭代次数矩阵,
Figure BDA0002822300630000082
表示矩阵b的奇异值,
Figure BDA0002822300630000083
表示矩阵A的奇异值
Figure BDA0002822300630000084
表示矩阵Ak,λ广义逆的奇异值,Ak,λ表示以k,λ为变量参数的矩阵。
如果迭代达到终止迭代次数k,GCV函数(即式15)达到最小值出现,或者达到残差的容差,此时则终止混合LSMR方法的迭代。
具体实施时,通过混合LSMR方法反演求解所述关系方程后得到海水表面的上行波波场,即可得到图2中B点的上行波波场U,再利用上行波拉东反变换算子Gu将所得的海水表面处上行波场进行延拓,进而得到检波端处波场,再进行FFT反变换,即可得到鬼波压制后的一次波结果Pz+Δz
Pz+Δz=Pze±iθ (17)
上式中:Pz表示海水表面处的波场,Δz表示z方向的采样间隔,
Figure BDA0002822300630000085
其中ω表示圆频率,v代表水速,ω=2πf,kx表示x方向的波数。
可见,上述鬼波衰减方法利用混合LSMR算法精确求解τ-p线性方程得到鬼波压制后的一次波波场,提高计算效率的同时改善了鬼波的压制效果;考虑鬼波延迟时间受出射角和偏移距的影响,弥补鬼波延迟时间估计中存在的误差,无需通过反演迭代求取最优鬼波延迟时间,提高了运算效率;鬼波衰减过程中人工干预减少,相对于目前现有鬼波衰减方法,具有较高的计算效率。
具体实施时,图3中(a)图是为两层水平层状模型,第一层速度为1500m/s,第二层速度为2000m/s。采用有限差分声波正演模拟,可以得到单炮带有鬼波信号的地震记录如图3中的(b)图所示,正演模拟时采用主频为350Hz的雷克子波,记录道数为101道,道间距为25m,采样间隔为1ms,记录长度为5s,单边放炮,拖缆深度为5~50m。
图4为图3中(b)图的单炮带有鬼波信号的地震记录进行压制鬼波后的单炮记录。其中,图4中(a)图显示的是采用传统LSMR算法压制鬼波后的单炮记录,图4中(b)图显示的是采用鬼波衰减方法的混合LSMR算法压制鬼波后的单炮记录,对比分析图4中的(a)图和(b)图可以看出,常规LSMR算法压制鬼波后的单炮记录鬼波仍有残余,上述鬼波衰减方法压制鬼波后的单炮记录明显优于常规LSMR算法。图5中(a)图为图3中单炮带有鬼波信号的地震记录数据对应的频谱图,图5中(b)图为图4中(a)图显示的采用传统LSMR算法压制鬼波后的单炮记录的频谱图,图5中(c)图为图4中(b)图显示的采用鬼波衰减方法的混合LSMR算法压制鬼波后的单炮记录的频谱图,可以看出,两种压制鬼波方法都能有效消除鬼波所引起的陷频特征,相比于常规LSMR方法,上述鬼波衰减方法得到的振幅谱中多个陷波频率的能量得到补偿,整个剖面高频、低频能量都得到提高。
具体实施时,图6和图7是实际数据测试效果,采用某深水区实际数据进行了应用测试。该工区为常规水平拖缆采集数据,气枪沉放深度7米,电缆深度8米,水深1000-2000米。采用上述鬼波衰减方法对地震数据进行鬼波压制。
图6测试了不同正则化参数对鬼波压制效果的影响,图6中(a)为鬼波压制前的地震数据,图6中(b)图为基于偏差原理(DP)方法确定正则化参数并应用于鬼波压制后得到的地震数据,图6中(c)图为基于无偏预测风险估计(UPRE)方法确定正则化参数并应用于鬼波压制后得到的地震数据,图6中(d)图为基于广义交叉验证(GCV)方法确定正则化参数并应用于鬼波压制后得到的地震数据。通过对比发现,广义交叉验证(GCV)方法相比偏差原理(DP)方法和无偏预测风险估计(UPRE)方法,对鬼波的压制效果更好,说明不同的正则化参数对鬼波压制起着重要的作用。
图7中的(a)图是鬼波压制前的单炮记录,图7中的(b)图是鬼波压制后的单炮记录。从图中对比可以看出,跟随在有效反射波之后鬼波(图7中(a)图内的箭头所指)得到了有效压制,反射波的相位由原来的‘白黑白’三个相位变成单一的‘白色相位’(图7中(b)图内箭头所指白色相位)。可见上述鬼波衰减方法对鬼波具有较好的压制作用,反射特征更加清晰,信噪比和分辨率都得到提升。
可见,上述鬼波衰减方法基于混合LSMR算法精确求解τ-p线性方程得到鬼波压制后的波场,提高计算效率的同时改善了鬼波的压制效果。通过利用合成数据对该鬼波衰减方法进行测试,验证了该鬼波衰减方法的有效性和适用性,对于提升海洋地震资料处理质量有着重要的潜在意义和实用价值。
在本实施例中,提供了一种计算机设备,如图8所示,包括存储器802、处理器804及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述任意的鬼波衰减方法。
具体的,该计算机设备可以是计算机终端、服务器或者类似的运算装置。
在本实施例中,提供了一种计算机可读存储介质,所述计算机可读存储介质存储有执行上述任意的鬼波衰减方法的计算机程序。
具体的,计算机可读存储介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。计算机可读存储介质的例子包括,但不限于相变内存(PRAM)、静态随机存取存储器(SRAM)、动态随机存取存储器(DRAM)、其他类型的随机存取存储器(RAM)、只读存储器(ROM)、电可擦除可编程只读存储器(EEPROM)、快闪记忆体或其他内存技术、只读光盘只读存储器(CD-ROM)、数字多功能光盘(DVD)或其他光学存储、磁盒式磁带,磁带磁磁盘存储或其他磁性存储设备或任何其他非传输介质,可用于存储可以被计算设备访问的信息。按照本文中的界定,计算机可读存储介质不包括暂存电脑可读媒体(transitory media),如调制的数据信号和载波。
基于同一发明构思,本发明实施例中还提供了一种鬼波衰减装置,如下面的实施例所述。由于鬼波衰减装置解决问题的原理与鬼波衰减方法相似,因此鬼波衰减装置的实施可以参见鬼波衰减方法的实施,重复之处不再赘述。以下所使用的,术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现,但是硬件,或者软件和硬件的组合的实现也是可能并被构想的。
图9是本发明实施例的鬼波衰减装置的一种结构框图,如图9所示,该装置包括:
海上地震资料获取模块902,用于获取海上水平拖缆地震资料;
方程建立模块904,用于基于平面波传播原理,根据所述海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
方程求解模块906,用于利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
数据处理模块908,用于对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
在一个实施例中,所述方程建立模块,包括:
第一正则化参数确定模块,用于基于最小化函数,通过以下公式确定正则化参数;
Figure BDA0002822300630000101
其中,GA,b(λ)表示正则化参数矩阵;λ表示正则化参数;n表示迭代总次数;b表示实测数据;A表示总波场的矩阵;I表示单位矩阵;
Figure BDA0002822300630000111
表示矩阵A的广义逆;xλ表示模拟数据;
Figure BDA0002822300630000112
表示矩阵A和模拟数据组成的矩阵。
在一个实施例中,还包括:
第二正则化参数确定模块,用于通过以下公式确定正则化参数:
Figure BDA0002822300630000113
其中,λ表示正则化参数;b表示实测数据;||e||表示噪声;xλ表示模拟数据;
Figure BDA0002822300630000114
表示矩阵A和模拟数据组成的矩阵。
在一个实施例中,还包括:
第三正则化参数确定模块,用于通过以下公式确定正则化参数:
Figure BDA0002822300630000115
其中,UA,b(λ)表示待求的正则化参数矩阵;λ为正则化参数;b表示实测数据;n表示迭代总次数;σ2表示噪声方差;A表示总波场的矩阵;
Figure BDA0002822300630000116
代表矩阵A的广义逆;xλ表示模拟数据;
Figure BDA0002822300630000117
表示矩阵A和模拟数据组成的矩阵。
在一个实施例中,还包括:
迭代终止模块,用于在最小平方残差方法迭代过程中,确定终止迭代次数,当迭代次数达到所述终止迭代次数时,终止迭代。
在一个实施例中,所述迭代终止模块通过以下公式确定终止迭代次数:
Figure BDA0002822300630000118
其中,G(k)表示终止迭代次数矩阵;k表示终止迭代次数;λ表示正则化参数;I表示单位矩阵;n表示迭代总次数;
Figure BDA0002822300630000119
表示矩阵b的奇异值;
Figure BDA00028223006300001110
表示矩阵A的奇异值;
Figure BDA00028223006300001111
表示矩阵Ak,λ广义逆的奇异值;Ak,λ表示以k,λ为变量参数的矩阵。
本发明实施例实现了如下技术效果:基于平面波传播原理,根据海上水平拖缆地震资料建立τ-p域拖缆观测的总波场与海水表面上行波波场之间的关系方程后,利用最小平方残差方法反演求解关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场,进而对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。提出了在最小平方残差方法迭代过程加入正则化参数,将正则化参数应用于最小平方残差方法的子问题,实现了采用混合最小平方残差方法来反演求解关系方程,有利于使得迭代解逼近真实解,进而有利于提高反演求解结果的准确性,有利于提高计算效率的同时改善鬼波的压制效果。
显然,本领域的技术人员应该明白,上述的本发明实施例的各模块或各步骤可以用通用的计算装置来实现,它们可以集中在单个的计算装置上,或者分布在多个计算装置所组成的网络上,可选地,它们可以用计算装置可执行的程序代码来实现,从而,可以将它们存储在存储装置中由计算装置来执行,并且在某些情况下,可以以不同于此处的顺序执行所示出或描述的步骤,或者将它们分别制作成各个集成电路模块,或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样,本发明实施例不限制于任何特定的硬件和软件结合。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (14)

1.一种鬼波衰减方法,其特征在于,包括:
获取海上水平拖缆地震资料;
基于平面波传播原理,根据所述海上水平拖缆地震资料建立拉东域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
2.如权利要求1所述的鬼波衰减方法,其特征在于,在每次迭代中确定正则化参数,包括:
基于最小化函数,通过以下公式确定正则化参数;
Figure FDA0002822300620000011
其中,GA,b(λ)表示正则化参数矩阵;λ表示正则化参数;n表示迭代总次数;b表示实测数据;A表示总波场的矩阵;I表示单位矩阵;
Figure FDA0002822300620000012
表示矩阵A的广义逆;xλ表示模拟数据;
Figure FDA0002822300620000013
表示矩阵A和模拟数据组成的矩阵。
3.如权利要求1所述的鬼波衰减方法,其特征在于,在每次迭代中确定正则化参数,包括:
通过以下公式确定正则化参数:
Figure FDA0002822300620000014
其中,λ表示正则化参数;b表示实测数据;||e||表示噪声;xλ表示模拟数据;
Figure FDA0002822300620000015
表示矩阵A和模拟数据组成的矩阵。
4.如权利要求1所述的鬼波衰减方法,其特征在于,在每次迭代中确定正则化参数,包括:
通过以下公式确定正则化参数:
Figure FDA0002822300620000021
其中,UA,b(λ)表示待求的正则化参数矩阵;λ为正则化参数;b表示实测数据;n表示迭代总次数;σ2表示噪声方差;A表示总波场的矩阵;
Figure FDA0002822300620000022
代表矩阵A的广义逆;xλ表示模拟数据;
Figure FDA0002822300620000023
表示矩阵A和模拟数据组成的矩阵。
5.如权利要求1至4中任一项所述的鬼波衰减方法,其特征在于,还包括:
在最小平方残差方法迭代过程中,确定终止迭代次数,当迭代次数达到所述终止迭代次数时,终止迭代。
6.如权利要求5所述的鬼波衰减方法,其特征在于,确定终止迭代次数,包括:
通过以下公式确定终止迭代次数:
Figure FDA0002822300620000024
其中,G(k)表示终止迭代次数矩阵;k表示终止迭代次数;λ表示正则化参数;I表示单位矩阵;n表示迭代总次数;
Figure FDA0002822300620000025
表示矩阵b的奇异值;
Figure FDA0002822300620000026
表示矩阵A的奇异值;
Figure FDA0002822300620000027
表示矩阵Ak,λ广义逆的奇异值;Ak,λ表示以k,λ为变量参数的矩阵。
7.一种鬼波衰减装置,其特征在于,包括:
海上地震资料获取模块,用于获取海上水平拖缆地震资料;
方程建立模块,用于基于平面波传播原理,根据所述海上水平拖缆地震资料建立拉东域拖缆观测的总波场与海水表面上行波波场之间的关系方程;
方程求解模块,用于利用最小平方残差方法反演求解所述关系方程,在最小平方残差方法迭代过程中,在每次迭代中确定正则化参数,并将正则化参数应用于最小平方残差方法的子问题,得到海水表面的上行波波场;
数据处理模块,用于对海水表面的上行波波场进行数据处理得到鬼波压制后的一次波。
8.如权利要求7所述的鬼波衰减装置,其特征在于,所述方程建立模块,包括:
第一正则化参数确定模块,用于基于最小化函数,通过以下公式确定正则化参数;
Figure FDA0002822300620000028
其中,GA,b(λ)表示正则化参数矩阵;λ表示正则化参数;n表示迭代总次数;b表示实测数据;A表示总波场的矩阵;I表示单位矩阵;
Figure FDA0002822300620000031
表示矩阵Aλ的广义逆;xλ表示模拟数据;
Figure FDA0002822300620000032
表示矩阵A和模拟数据组成的矩阵。
9.如权利要求7所述的鬼波衰减装置,其特征在于,所述方程建立模块,还包括:
第二正则化参数确定模块,用于通过以下公式确定正则化参数:
Figure FDA0002822300620000033
其中,λ表示正则化参数;b表示实测数据;||e||表示噪声;xλ表示模拟数据;
Figure FDA0002822300620000034
表示矩阵A和模拟数据组成的矩阵。
10.如权利要求7所述的鬼波衰减装置,其特征在于,所述方程建立模块,还包括:
第三正则化参数确定模块,用于通过以下公式确定正则化参数:
Figure FDA0002822300620000035
其中,UA,b(λ)表示待求的正则化参数矩阵;λ为正则化参数;b表示实测数据;n表示迭代总次数;σ2表示噪声方差;A表示总波场的矩阵;
Figure FDA0002822300620000036
代表矩阵A的广义逆;xλ表示模拟数据;
Figure FDA0002822300620000037
表示矩阵A和模拟数据组成的矩阵。
11.如权利要求7至10中任一项所述的鬼波衰减装置,其特征在于,还包括:
迭代终止模块,用于在最小平方残差方法迭代过程中,确定终止迭代次数,当迭代次数达到所述终止迭代次数时,终止迭代。
12.如权利要求11所述的鬼波衰减装置,其特征在于,所述迭代终止模块通过以下公式确定终止迭代次数:
Figure FDA0002822300620000038
其中,G(k)表示终止迭代次数矩阵;k表示终止迭代次数;λ表示正则化参数;I表示单位矩阵;n表示迭代总次数;
Figure FDA0002822300620000039
表示矩阵b的奇异值;
Figure FDA00028223006200000310
表示矩阵A的奇异值;
Figure FDA00028223006200000311
表示矩阵Ak,λ广义逆的奇异值;Ak,λ表示以k,λ为变量参数的矩阵。
13.一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现权利要求1至6中任一项所述的鬼波衰减方法。
14.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有执行权利要求1至6中任一项所述的鬼波衰减方法的计算机程序。
CN202011420878.2A 2020-12-08 2020-12-08 鬼波衰减方法、装置、计算机设备及可读存储介质 Pending CN114609671A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011420878.2A CN114609671A (zh) 2020-12-08 2020-12-08 鬼波衰减方法、装置、计算机设备及可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011420878.2A CN114609671A (zh) 2020-12-08 2020-12-08 鬼波衰减方法、装置、计算机设备及可读存储介质

Publications (1)

Publication Number Publication Date
CN114609671A true CN114609671A (zh) 2022-06-10

Family

ID=81857033

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011420878.2A Pending CN114609671A (zh) 2020-12-08 2020-12-08 鬼波衰减方法、装置、计算机设备及可读存储介质

Country Status (1)

Country Link
CN (1) CN114609671A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117148443A (zh) * 2023-10-27 2023-12-01 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117148443A (zh) * 2023-10-27 2023-12-01 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法
CN117148443B (zh) * 2023-10-27 2024-03-19 胜利信科(山东)勘察测绘有限公司 基于鬼波提取与转换的浅剖数据信噪比增强方法

Similar Documents

Publication Publication Date Title
US9625593B2 (en) Seismic data processing
US7599798B2 (en) Migrating composite seismic response data to produce a representation of a seismic volume
US20120275267A1 (en) Seismic Data Processing
US8902705B2 (en) Regularisation of irregularly sampled seismic data
US10169652B2 (en) Spatial expansion seismic data processing method and apparatus
Huang et al. Source-independent extended waveform inversion based on space-time source extension: Frequency-domain implementation
US10215869B2 (en) System and method of estimating anisotropy properties of geological formations using a self-adjoint pseudoacoustic wave propagator
US20170031045A1 (en) Method and apparatus for modeling and separation of primaries and multiples using multi-order green's function
GB2547942A (en) Method for deghosting and redatuming operator estimation
CN109946741B (zh) 一种TTI介质中纯qP波最小二乘逆时偏移成像方法
US10795039B2 (en) Generating pseudo pressure wavefields utilizing a warping attribute
Majdański et al. Attenuation of free-surface multiples by up/down deconvolution for marine towed-streamer data
Yilmaz et al. Discrete plane-wave decomposition by least-mean-square-error method
Chemingui et al. Seismic data reconstruction by inversion to common offset
CN111665556B (zh) 地层声波传播速度模型构建方法
AU2008327686A1 (en) Moveout correction of seismic data
CN114609671A (zh) 鬼波衰减方法、装置、计算机设备及可读存储介质
Operto et al. 3D ray+ Born migration/inversion—Part 2: Application to the SEG/EAGE overthrust experiment
CN117388920A (zh) 一种基于迭代地震干涉法的近偏移距数据重建方法
NO348176B1 (en) Improvement to seismic processing based on predictive deconvolution
de Jonge et al. Deghosting dual‐component streamer data using demigration‐based supervised learning
CN109425892B (zh) 地震子波的估计方法及系统
CN111538088B (zh) 一种海上斜缆波场校正方法
CN115061197A (zh) 二维海面鬼波水体成像测量方法、系统、终端及测流设备
Ernst et al. Reduction of near-surface scattering effects in seismic data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination