CN113945982A - 用于去除低频和低波数噪声以生成增强图像的方法和系统 - Google Patents

用于去除低频和低波数噪声以生成增强图像的方法和系统 Download PDF

Info

Publication number
CN113945982A
CN113945982A CN202110734440.XA CN202110734440A CN113945982A CN 113945982 A CN113945982 A CN 113945982A CN 202110734440 A CN202110734440 A CN 202110734440A CN 113945982 A CN113945982 A CN 113945982A
Authority
CN
China
Prior art keywords
transitory computer
readable device
domain
retrieved
survey area
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.)
Granted
Application number
CN202110734440.XA
Other languages
English (en)
Other versions
CN113945982B (zh
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.)
China Petroleum and Chemical Corp
Original Assignee
China Petroleum and Chemical Corp
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 China Petroleum and Chemical Corp filed Critical China Petroleum and Chemical Corp
Publication of CN113945982A publication Critical patent/CN113945982A/zh
Application granted granted Critical
Publication of CN113945982B publication Critical patent/CN113945982B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G01V1/50—Analysing data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6169—Data from specific type of measurement using well-logging
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/70—Other details related to processing
    • G01V2210/74—Visualisation of seismic data

Abstract

公开了一种计算机实现的逆时偏移方法，其中使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，从而生成增强的图像。

Description

用于去除低频和低波数噪声以生成增强图像的方法和系统
技术领域
本发明大体上涉及通过将输入地震道数据经偏移算子映射到输出结构图像来对复杂的地下结构进行成像，特别是使用3D拉普拉斯算子并结合2D滤波技术，以去除低频和低波数噪声，从而阐明和解释地下碳酸氢盐目标以用于地震后处理。
背景技术
1.引言
本领域的普通技术人员已经知道，岩性识别和流体鉴定在地震勘探和储层表征中起着重要的作用。岩性识别和流体鉴定的定义是使用叠前地震数据的加权叠加方法来估算的，最早由Smith,G.和Gidlow,P.在他们的研究“Weighted stacking for rock propertyestimation and detection of gas:Geophysical Prospecting(1987)”中提出。后来，Goodway，B.、T.Chen和J.Downton在SEG第67届国际年会上介绍了“Improved AVO fluiddetection and lithology discrimination using Lamé petrophysical parameters:“λρ”,“μρ”&“λ/μfluid stack”from P and S inversions”，发现拉梅参数对烃饱和度来说比弹性参数更敏感，并提出了用于岩性识别和流体鉴定的LMR方法。之后发现，该方法会取决于埋藏深度、压实度和上覆压力而发生变化，因此对于具有多孔介质的实际储层而言有一定的限制。在2006年，Quakenbush,M.、B.Shang和C.Tuttle引入了泊松阻抗的概念，将其用于岩性识别和流体鉴定。实际上，在他们的被认为较新颖的研究论文“Poisson impedance:The Leading Edge,v.25,no.2,pp.128-138”中，首次根据区分岩性和流体含量的所述方法开发了一种敏感性工具。这开启了开发计算机实现的方法以获取流体属性、执行流体鉴定但不包括烃预测的时代。因此，在1980年代或前后，逆时偏移迅速成为复杂地质学中首选的高端成像方法(例如参见Baysal,E.、D.D.Kosloff和J.W.C.Sherwood,1983,Reverse timemigration:Geophysics,Vol.48,P1514-1524；以及Etgen,J.,1986,Prestack reversetime migration of shot profiles:SEP-Report,Vol.50,P151-170)。
然而，RTM的计算挑战之一是震源波场在时间上向前传播(从0到t时间最大值)，而接收器波场在时间上向后传播(从t时间最大值到0)，但成像步骤要求这两个波场在每个时间t处均相关联。人们很快发现，将波场之一存储在中央处理单元或图形处理单元(内存资源)中对于3D问题来说是不切实际的。最明显的解决方案是将每个成像步骤的波场存储在物理磁盘上，例如硬盘。这又需要每个节点上的大量磁盘存储，并导致问题迅速成为输入/输出(IO)限制。为了解决这个问题，提出了一种检查点方案，其中将较少量的快照存储到硬盘驱动器中，并重新生成中间波场(参见Symes,W.W.,2007,Reverse time migration withbest checkpointing:Geophysics,Vol.72,P SM213-SM221)。另一种提出的方法是保存计算域边界处的波场并重新注入它们(参见Dussaud,E.,W.W.Symes,L.Lemaistre,P.Singer,B.Denel以及A.Cherrett,2008,Computational strategies for reverse-timemigration:78th Annual Internat.Mtg.,Soc.Expl.Geophys.,Expanded Abstracts,PSPMI 3.3)。
另一方面，根据地震振幅以及振幅-偏移关系来进行烃预测仍然是一项艰巨的任务。一种方法是使用地震反射来将其与地下岩石特性紧密相关。然而，地震数据中最强的振幅-偏移关系(AVO)通常是由岩石中的烃饱和引起的。使用叠前地震反演来提取地震数据的关于地下弹性参数方面的信息的进步极大地帮助了以最小的误差来表征岩相和预测储层属性，从而减少了世界上某些盆地中的干井数量和钻井风险(例如参见Russel,B.，2014，Prestack seismic amplitude analysis:an integrated overview:Interpretation,v.2,no.2,SC19-SC36)。这种叠前地震反演模型已经常规地应用于岩性预测和流体探测，以识别油气勘探中的潜在目标。最近，它已被广泛用于估算非常规页岩气应用中的甜点。
叠前地震反演模型的本质在于这样的事实：当流体饱和度改变时，岩石的剪切模量不会改变。当流体饱和度改变时，体积模量则发生显著变化。由于岩石骨架的剪切模量等于具有孔隙流体的岩石的剪切模量，因此可以将剪切阻抗视为主要连接到岩石骨架的地震属性，而声阻抗则由孔隙流体和岩石骨架两者主导。结果，泊松阻抗可以最佳地从剪切模量和声阻抗中消除岩石骨架的影响，从而可以产生更好的流体含量的分辨率。通过叠前地震数据，传统上是根据可以直接从地震数据中求反的P波速度、S波速度和密度而间接地计算出泊松阻抗(见Goodway，同上)。但是，参数估算的间接方法会因间接计算而产生更多的不确定性(例如参见Zong,Z.、X.Yin和G.Wu,2013,Elastic impedance parameterizationand inversion with Young’s modulus and Poisson’s ratio，Geophysics，v.78,no.6,N35-N42)。在2006年，Wang,B.、X.Yin和F.Zhang在他们的研究论文“Lamé parametersinversion based on elastic impedance and its application:Applied Geophysics,v.3,pp.120-123”中将从叠前地震数据中估计拉梅参数的直接方法与使用反演P波和S波来估计拉梅参数的间接方法进行了比较，但是得出的结论是，它仍然会造成更大的偏差。
因此，由于所有这些试图减少累积误差以提供准确而合理的结果的尝试均告失败，直接和间接地从叠前地震反演中提取流体因子近年来受到了广泛的关注，特别是考虑到解决上述大型科学和工程问题所需要的具有高性能计算系统的系统计算机的更广泛的可用性和可得性。
2.使用拉普拉斯算子的逆时偏移
如前所述，逆时偏移(RTM)是一种先进的地震偏移成像方法。然而，当震源波场和接收器波场互相关时，直接到达波、背向散射波和翻转波的互相关会产生低频和低波数的噪声，从而掩盖了最终的成像结果。因此，在文献中已经提出了几种技术，例如速度平滑高通滤波器(参见Mulder,W.A.和Plessix,R.-E.,2003,One-way and two-way waveequation migration,Vol.73,SEG Annual Meeting,P881-884)、Poyntingvectors(参见Yoon,K.、Marfurt,KJ和Starr,W.,2004,Challenges in reverse-time migration,Vol.74,SEG Annual Meeting,P1057-1060)，以及界面上的方向阻尼终端(参见Fletcher,RF、Fowler,P.、Kitchenside,P.和Albertin,U.,2005,Suppressing artifacts inprestack reverse time migration,Vol.73,SEG Annual Meeting,P2049-2051)。在实践中，这些方法难以正确实施，或者具有会不期望地扭曲偏移图像的频谱或幅度的显著缺点。后来又提出发一种新的成像条件(见Liu,F.、Zhang,G.、Morton,S和Leveille,J.，2007，Reverse-time migration using oneway wavefield imaging conditions,Vol.77,SEGAnnual Meeting,P2170-2174)，以解决现有技术提到的问题。它通常包括首先将波场分解为单向分量，然后仅对作为反射出现的波分量进行互相关。然而，这种技术在2D空间域中最有效，因为在3D空间域中完全分解波场是计算密集型的，并且还包含有不完全的分解，导致无意中去除了复杂结构中的陡峭倾斜的反射体。
因此，引入3D空间域中的拉普拉斯滤波，作为抑制这些低频和低波数噪声的常用方法，以适应任何复杂的介质，但这会增加计算成本。然而，后来发现简单、直接的拉普拉斯滤波会破坏从勘探区域获取的有用信号的特征。Zhang和Sun(参见Yu Zhang和James Sun,January 2009,Practical issues in reverse time migration:true amplitudegathers,noise removal and harmonic source encoding,First Break,V26,P29-37)提出了一种改进的拉普拉斯滤波方法，并给出了简单的计算公式和说明。该方法在保留有用信号特征的同时有效地抑制了逆时偏移中的低频噪声，但缺乏详细和严格的数学推导。
后来，引入了几种基于带有分数拉普拉斯算子的粘声波动方程的逆时偏移(RTM)方案来补偿粘性效应。在这些方法中，首先发展了一种基于粘声波动方程的高效的波传播模拟方法；然后将粘声波动方程分解为幅度衰减方程和相位分散方程，以实现更合理的Q补偿RTM算法；最后使用简单的地堑模型和更逼真的修正Marmousi模型来测试Q补偿的逆时偏移方法。然而，这些方法提供的补偿偏移结果仍然非常接近于没有衰减的地震数据的常规RTM获得的结果。
因此，针对现有技术的一些缺陷，本发明实施例引入了一种新的计算机实现方法，用于去除低频和低波数噪声，以生成增强图像。
3.低通波数滤波(k-filtering)
k滤波技术是一种根据频率和波向量中的波或冲击能量分布来表征平稳波动的方法。简单地说，它可以被认为是定位特定信号的速度和方向。时空信号滤波器的设计遵循与一维信号类似的方法。一维信号滤波器设计成：如果滤波器的要求是提取特定的非零频率范围内的频率分量，则可以确定具有适当通带和阻带频率的带通滤波器。类似地，在多维系统的情况下，滤波器的波数-频率响应设计为在(k,ω)a.k.a.波数-频率和其他地方为零的设计区域内统一。因此，低通波数滤波具有区分相同频率的波模式，甚至具有不同波矢的波模式的能力。该方法基于波场分量的同时多点测量，其中使用滤波器组来提高分辨率。k滤波技术最突出的用途以前仅用于磁场测量。
事实上，k滤波技术是一种多航天器方法，用于估计波场能量分布P，它是角频率ω和波矢量k的函数。k滤波的基本思想是将空间域中的滤波器组应用于在多个空间位置测量并已转换到频域的数据时间序列(参见Pincon,J.-L.和U.Motschmann(1998)，Generalframework,in Analysis Methods for Multi-Spacecraft Data,edited by G.Paschmannand P.W.Daly,ISSI Sci.Rep.,SR-001,chap.3,pp.65–78)。这个特定的滤波器取决于ω和k，其理想设计方式是吸收所有波场能量，除了对应于具有ω和k的平面波的部分，它应该不受干扰地通过滤波器。
k滤波技术可以被认为是一种将总测量空间相关矩阵M(ω)分解为相关矩阵之和的方法，对应于纯平面波模式。
M(ω)＝[A(ω)A'(ω)] (1)
其中，
Figure BDA0003139878800000051
在空间点r处。
由于M(w)是一个NL×NL矩阵，因此它只能分解为NL个线性无关的相关矩阵。这些矩阵之一对应于数据中的非相干噪声，这使得可以使用k滤波方法区分理论最大值为NL-1的不同种波模式。通过测量波场的更多分量，本领域的普通技术人员可以区分更多种的波模式。使用更多数量的接收传感器也会有所帮助，因为波模式的数量取决于测量空间相关矩阵的大小。另一方面，使用场分量的约束将减少波模式的数量，因为它们减少了波场的独立分量的数量。因此，k滤波方法可以通过方便地选择约束矩阵来抑制波场能量分布中的任何混叠峰值。当测量更多的波场分量时，这可以更有效地完成，但它仍然与执行某些抗混叠有关。
4.抗混叠
地震偏移的过程包括通过偏移算子将输入地震道数据映射到输出结构图像。偏移过程的空间混叠在所有三个域中都可能存在问题：数据、算子和图像。这三种类型的空间混叠是不同且独立的，尽管它们经常被混淆。
首先，当记录的地震数据中沿反射和衍射事件存在违反了以下标准的地震频率f和倾角θ时，会发生数据混叠：
Figure BDA0003139878800000052
参数fd是针对给定轨迹间距Δx、接收器位置处的速度vr和记录表面处的局部平面波入射角θr的数据中的最大非混叠频率。一旦数据在记录过程中出现混叠，就很难避免后续的处理伪影，除非以更精细的组间距重新勘探，或者尝试在不插入伪影的情况下插入原始数据(例如参见Spitz,S.,1991,Seismic trace interpolation in the F-X domain:Geophysics,V56,P785-794；或者Claerbout,J.F.,1992,Earth soundings analysis:Processing versus inversion:Blackwell Scientific Publications)。
另一方面，当图像空间的输出采样太粗糙而无法正确表示偏移的倾角时，就会出现图像混叠。在基尔霍夫偏移的情况下，图像混叠主要是一个表面问题，可以通过选择合适的图像网格来解决，以估计偏移的结构图像。请注意，混叠图像空间中每个网格点的基尔霍夫偏移值是未受污染且正确的(假设没有数据或运算符混叠)。例如，井位置处的单个偏移输出迹线可能是完全准确的，而无需估计任何相邻的偏移迹线。这种基尔霍夫偏移的“面向目标”能力是其他类型的偏移方法(如光谱或有限差分)无法实现的。因此，如果需要将偏移图像作为后续多通道处理的输入，图像混叠是一个问题。
然后是算子混叠，当算子沿偏移求和轨迹的倾斜对于给定的输入地震道间距和频率内容来说太陡峭时，就会发生这种情况。基尔霍夫积分算子的空间混叠很可能是因为这些算子是在离散时空域中设计的，而不是在谱域中设计的。必须非常小心以确保给定时空算子的f-k谱不包含超过奈奎斯特频率和波数的混叠能量。这意味着空间欠采样并包含显著陡倾和高频能量的时空算子可能会遭受严重的算子混叠。相比之下，在f-k域中明确设计的偏移算子根据定义是无混叠的，而在f-k域中设计的算子很容易修改，以抑制空间混叠。
通常，抗混叠算子可以通过轨迹插值、孔径控制和算子倾角滤波来部分地实现。迹线插值减少了空间采样间隔，从而增加了空间奈奎斯特频率，但增加了总3-D偏移输入负载的潜在成本。偏移孔径控制和算子倾角滤波具有抑制对输出图像的混叠陡倾角贡献的效果。不幸的是，在对盐翼和断层等近垂直结构进行成像时，这些被抑制的陡倾可能正是感兴趣的事件。之后提出了一种抗混叠倾角动校正(DMO)算子的方法(SeeHale,D.,1991,Anonaliasedintegral method for dip moveout:Geophysics,V56,P795-805)，其中复制输入地震道的时移版本，以沿着算子的等时采样版本来近似求和。这种空间插值意味着沿时间方向对输入轨迹值进行boxcar平均，以确保抑制算子混叠和图像混叠，这会调用比对于单独偏移算子抗混叠的情况来说是绝对必要的更多的低通滤波和分辨率损失。因此，从浮点运算和内存的角度来看，这种时移输入轨迹和空间插值的复制被证明对于同时处理数千条轨迹的3D偏移算法来说是低效的。
为了解决上文中的Hale D.冲突，随后提出了一种针对抗混叠偏移算子的近似但具有成本效益的方法(参见Gray,SH,1992，Frequency-selective design of theKirchhoff migration operator:Geophysical Prospecting,V40,P565-571)。Gray S.H.,提议通过在不同的高频截止频率下对输入数据的大约三个副本进行低通滤波来对输入轨迹数据执行依赖于孔径的低通时间滤波。在偏移期间，偏移算子的大孔径部分将越来越多地从最低频率内容过滤的轨迹副本中提取。这对应于这样的概念：偏移算子的大孔径部分具有最陡峭的倾角，因此只有最低频率将不混叠。
无论如何，算子抗混叠标准通常被设置为在算子轨迹与地震道相交的位置处定义地震数据中的最大非混叠时间频率，并且是局部算子倾角和道间距的函数。为了保持无混叠，沿着算子轨迹求和的地震道样本的序列必须满足以下倾角的奈奎斯特采样标准：
Figure BDA0003139878800000071
因子ΔT是沿相邻轨迹之间的偏移算子的时差时间。术语
Figure BDA0003139878800000072
是偏移算子轨迹在与地震道相交点处的局部空间导数，而Δρ是沿记录面的地震道间距。抗混叠方程(2)建议对地震道进行局部低通滤波，以确保大于fmax的频率不会错误地偏移到输出图像中。
尽管如此，这些现有的滤波方法主要用于时间到深度(或相反)的转换，它们发生在偏移过程中，它们都会导致图像质量差和模糊。因此，所提出的计算机实现方法实现了已被证明更准确的方法的组合，因为它们不仅为每个衍射轨迹分别计算高截止频率的完整向量，而且因为它具有以下成本效益：不需要大量数据副本的内存密集型存储来实现抗混叠。
5.层析反演
深度偏移道集的层析成像是一种改进速度-深度模型的方法。层析成像的原理是，如果使用正确的速度深度模型进行偏移，图像道集应该在所有接收传感器处具有相同的事件深度(Tian-wen Lo和Philip Inder Weisen，1994，Fundamentals of seismicTomography,SEG monograph series)。因此，层析成像试图通过分析PSDM后的残余延迟来纠正速度深度模型中的错误。当使用基于非全局方法的反演方法得出的初始错误速度模型来执行叠前深度偏移时，深度道集将表现出非平坦性。非平坦度是模型误差的度量。层析成像使用这种非平坦度的测量(剩余时差分析)作为输入，并试图找到一个可最大限度地减少误差的替代模型。层析成像原理将时间误差归因于速度和深度这两者的误差。
然而，在反演理论的领域中，对象是基于与该对象相关联的测量或观察来描述的。对于工业规模的应用，没有足够的数据来确定唯一的解决方案，并且所获得的数据可能有噪声和/或不可靠。因此，数学的整个分支都随着尝试根据对不准确、不充分和不一致的数据的解释来估计解决方案而发展。在地震实验中进行传播时间测量的情况下，试图解决的反演问题是地球的速度结构。
当将层析成像的概念与多维结构的反演合并时，可以观察到它自1970年代后期以来一直用于油田勘探。然而，这种反演的分辨率通常以公里甚至几十公里为单位。由于模型中包含了数千个节点，因此无法解决与油田勘探的复杂系统有关的细节，因此3D成像仍然基于集中式离线处理，并通常由多个有源记录完成，其中超过多年跨度的变化是主要目标。为了实现这一点，需要新的方案和方法来解决实时地震层析成像反演问题，这是使用单一计算算法无法实现的。这是开发可以传递准确成像系统和模型的计算机实现方法的一个动机。
由于层析成像基于射线追踪，因此可以针对反射、透射和折射进行公式化。在地震反射勘探中用于计算静力学校正的几种技术使用了折射层析成像，而透射层析成像用于跨井应用，其中震源和接收器都在介质内(在钻孔内)；因此，可以访问任何所传输的到达信息。除了到达时间之外，利用振幅信息可以进一步帮助基于射线的层析成像，以估计可靠的速度模型(例如参见Semtchenok NM、Popov MM和Verdel AR，2009，Gaussian beamtomography,Extended Abstracts,71st EAGE Conference&Exhibition,Amsterdam,U032)。除了速度估计之外，层析成像还可用于估计其他地球参数，例如吸收。此外，层析成像尝试使用穿过模型中单元的许多射线路径来求解一组联立方程。对于共中点(CMP)道集，给定的地下反射器元件有许多传播时间测量值，因此这五个射线路径的传播时间表达式可以写为：
T＝D*S (4)
其中，T是沿射线路径的总传播时间，D是速度单元中的射线路径的长度，S是速度单元中的慢度(即速度的倒数)。
不幸的是，在大多数情况下，矩阵D是不可逆的。要成为可逆的，它需要是方形的(即，传播时间测量的数量必须恰好与模型中的速度单元的数量相同)，并且还满足其他一些标准。因此，取而代之的是，该算法然后被反演为时间延迟和慢度更新，从而得到表达式：
ΔT＝D*ΔS (5)
然而，在该方法的上述描述中存在一些循环论证，因为要通过射线追踪来估计单元中的射线路径段长度，需要每个单元中反射器段的慢度更新和局部倾角估计。因此，当方法以正向建模开始时，通过使用模型的初始猜测的光线追踪，以及计算的相关联的传播时间，本发明提供第一猜测或估计模型。然后必须迭代层析成像，以收敛于真实模型的最佳估计，这通过最小化观察到的传播时间与通过射线追踪计算出的用于模型当前猜测的传播时间之间的差异来实现。话虽如此，本领域技术人员会意识到，算法(2)还有其他求解方式，例如：直接求解器，这是一种一步式解决方案，但仅适用于较小规模的问题；或迭代求解器，例如共轭梯度法(Scales,J.A.,1987，Tomographic inversion via the conjugategradient method,Geophysics,V52,P179-85)，它适用于大规模系统。
迄今为止描述的层析反演过程是根据未偏移数据的观察的和射线追踪的传播时间来进行的。然而，在过去十年中，本领域的技术人员主要处理深度偏移域中的数据和测量，因此还需要修改层析成像算法以考虑测量域的变化。因此，偏移域中的层析反演可以在叠前时间偏移之后进行(参见示例Pierre Hardy和Jean-Paul Jeannot，1999，3-Dreflection tomography in time-migrated space,SEG Technical Program ExpandedAbstracts,P1287-1290)或叠前深度偏移(例如Christof Stork,1992,Reflectiontomography in the post-migrated domain,Geophysics,V57,P680-692。)，在这两种情况下都产生深度模型，但最广泛的工业实践是使用叠前深度偏移数据的测量值进行反演。同样，当这些模型中的任何一个单独使用时，对于本领域的技术人员来说，它们对岩性识别的用处不大。
6.全波形反演
全波形反演(FWI)是非线性数据拟合程序，旨在从地震数据中获得地下特性的详细估计，这可以是被动或主动地震实验的结果。给定对地下参数的初始猜测(也称为模型)，通过求解波动方程来预测速度模型和各向异性参数模型。然后更新这些模型，以减少观察数据和预测数据之间的失配。这以迭代方式重复，直到数据失配足够小。
在数学上，FWI可以表述为M实验的优化问题：
Figure BDA0003139878800000091
其中m代表模型，
Figure BDA0003139878800000101
代表相位角，ρ(·)是惩罚函数。在惩罚中，di＝d1,…,dM是观测数据，Fii＝F1,…,FM是第i个实验的对应建模算子。优化问题通常使用以下形式的迭代算法来解决：
mk+1＝mk+αk×sk (7)
其中，sk称为搜索方向，αk是步长。最简单的例子是最速下降算法，在这种情况下，搜索方向由目标的负梯度给出：
Figure BDA0003139878800000102
如本领域中所观察到的那样，所有类型的波都参与到波的优化中，包括潜水波、超临界反射和多重散射波，例如多次波。用于正向建模的技术各不相同，包括：体积法，例如有限元法(参见Marfurt,K.,1984,Accuracy of finite-difference and finite-elementsmodeling of the scalar and elastic wave equation:Geophysics,V49,P533-549；Min,D.-J.、C.Shin、R.G.Pratt和H.S.Yoo,2003,Weighted-averaging finite-element methodfor 2D elastic wave equations in the frequency domain:Bulletin of theSeismological Society of America,V93,P904–921)、有限差分法(参见Virieux,J.,1986,P-SV wave propagation in heterogeneous media,velocitystress finitedifference method:Geophysics,V51,P889-901)、有限体积方法(Brossier,R.、J.Virieux和S.Operto,2008,Parsimonious finite-volume frequency-domain method for2D P-SVwave modelling:Geophysical Journal International,V175,P541-559)和伪谱方法(Danecek,P.和G.Seriani,2008,An efficient parallel Chebyshev pseudo-spectralmethod for large-scale 3D seismic forward modelling:70th Conference&TechnicalExhibition,EAGE,Extended Abstracts,P046)；边界积分法，例如反射率方法(Kennett,B.L.N.,1983,Seismic wave spread in layerified media:Cambridge UniversityPress)；广义屏幕法(Wu,R.S.,2003,Wave propagation,scattering and imaging usingdual-domain one-way and one-return propagation:Pure and Applied Geophysics,V160,P509-539)；离散波数法(Bouchon,M.、M.Campillo和S.Gaffet,1989,A boundaryintegral equation-Discrete wavenumber representation method to study wavepropagation in multilayered media having irregular interfaces:Geophysics,V54,P1134-1140)；广义射线法，例如WKBJ和Maslov地震图(Chapman,C.,1985,Ray theory andits extension:WKBJ and Maslov seismograms:Journal of Geophysics,V58,P27-43)；全波理论(de Hoop,A.,1960,A modification of Cagniard’s method for solvingseismic pulse problems:Applied Scientific Research,V8,P349-356)；以及衍射理论(Klem-Musatov,K.D.和A.M.Aizenberg,1985,Seismic modelling by methods of thetheory of edge waves:Journal of Geophysics,V57,P90-105)。然而，FWI尚未被认为是一种有效的地震成像技术，因为开创性的应用仅限于地震反射数据。对于短偏移距采集来说，地震波场对中间波长相当不敏感；因此，优化不能通过迭代更新来充分地重建真实的速度结构。只有当提供足够精确的初始模型时，波形拟合才能通过此类更新收敛到速度结构。为了对初始模型进行采样，已经使用全局和半全局技术进行了复杂的研究(例如参见Koren,Z.、K.Mosegaard、E.Landa、P.Thore和A.Tarantola,1991,Monte Carlo estimationand resolution analysis of seismic background velocities:Journal ofGeophysical Research,V96,P20289-20299；或者Mosegaard,K.和A.Tarantola,1995,Monte Carlo sampling of solutions to inverse problems:Journal of GeophysicalResearch,V100,P12431-12447)。由于对中间波长不敏感而导致的这些研究的性能相当差，这使许多研究人员相信这种优化技术并不是特别有效。
只有利用长偏移距和传输数据来重建结构的大中波长，FWI才算成熟(参见Mora,PR,1987,Nonlinear two-dimensional elastic inversion of multi-offset seismicdata:Geophysics,V52,P1211-1228；或者Pratt RG,1999,Seismic waveform inversionin the frequency domain,Part I:Theory and verification in a physical scalemodel:Geophysics,V64,P888-901页)。FWI试图在模型的每个点处表征一个广泛而连续的波数谱，将宏观模型构建和偏移任务重新统一到一个程序中。历史跨孔和广角地表数据示例说明了同时重建整个空间光谱的能力。然而，由于增加了在几十个波长和各种入射角上传播的波场所引入非线性性，FWI对长偏移距数据的稳健应用仍然具有挑战性(Sirgue L.,2006,The importance of low frequency and large offset in waveform inversion:68th Conference&Technical Exhibition,EAGE,Extended Abstracts,A037)。
发明内容
本发明公开了一种用于定位和成像地下反射体的新方法，其中在不同空间域中使用多步式计算机实现方法，以便通过平衡横向和垂直幅度来抑制低波数和低频噪声。该方法将产生位于勘探区域内的地下反射体的图像，其具有更高的横向分辨率和波数，以及在复杂介质中更高的高频截止频率和更低的低频截止频率，这是本领域的其他通常已知的方法无法实现的。因此，本发明的一个提议实施例包括在执行3D拉普拉斯运算之后仅在横向(x,y)方向上(即不在深度z方向上)执行低通k滤波。对于深度z方向，本发明的一个提议实施例使用提供的速度模型将其转换到时域，然后以高截止频率fmax执行低通频率滤波(例如f滤波)，其中所述fmax的值是用户预定义的输入参数。这样，在利用了速度模式的信息后，对于所采集的地震图像的所有深度来说，图像的频率内容都可以准确地包含在fmax的范围内并保持完整。此后将其转换回深度域，并生成勘探区域的3D空间域方向x、y和z上的最终堆叠图像。与任何现有方法相比，在复杂介质中，所述最终堆叠图像在横向分辨率和波数以及振幅平衡方面得到了增强。
下文将参考下面列出的附图描述本发明的进一步细节、示例和方面。
附图说明
通过结合附图并考虑以下描述，可以容易地理解本发明的教导。
图1示意性显示了根据本发明的某些实施例的说明性环境的横截面图的示意图，其包含地震源、地震接收器、井眼位置、井筒、各种传输射线以及各种入射角；
图2a和图2b是根据本发明的某些实施例的计算机实现的逆时偏移方法的流程图，其使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，以便生成增强的图像；
图3是数字高性能计算系统的方框形式的电气图，该数字高性能计算系统编程用于执行根据本发明的某些实施例的计算机实现的逆时偏移方法，其使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，以便生成增强的图像；和
图4示出了根据本发明的某些实施例的计算机实现的逆时偏移方法的图形用户界面。
具体实施方式
现在将详细参考本发明的几个实施例，其示例在附图中示出。应注意的是，在适用之处，在附图中使用相同或类似的附图标记来指示相同或类似的功能。附图仅出于说明的目的描绘了本发明的实施例。本领域的技术人员将从以下描述中容易地认识到，在不脱离本文描述的原理的情况下，可以采用其中示出的结构、系统和方法的替代实施例。
由于传统的逆时偏移方法通常直接在深度z方向上执行低通k滤波，并且从不使用任何关于速度变化的信息(它们隐含地假设所有三个空间方向(x,y,z)上的所有速度都是恒定的)，这些方法都不能将图像的频率内容包含在用户所需的频率变量fmax的精确范围内，从浅到深的任意位置上都是如此。因此，这些方法需要用户进行无休止的“试验和观察”经验计算，以得出z方向上适当的k滤波高截止参数。最终，这可能导致即使是本领域的普通技术人员也可能最终会损坏图像的频谱。因此，本文提出了创新的计算机实现方法的相关性。
图1是可以使用本发明的优选实施例的勘探区域101上方的地球的一部分的横截面图。重要的是应当注意，图1的勘探区域101是基于陆地的区域102。本领域的普通技术人员将认识到，地震勘探区域产生局部地质的详细图像，以确定可能的碳氢化合物(石油和天然气)储层的位置和大小，以及井位103。在这些勘探区域中，声波在爆炸期间从地下岩层在不同的入射点104处反弹，反射回地表的波由地震数据记录传感器105捕获，由数据传输系统305从所述传感器105处无线传输303，然后存储到内存资源304中以供后续处理，并且由非暂时性计算机可读设备306来处理，该计算机可读设备306由图3所示的计算机系统设备307来控制。
参照图1，勘探区域上方的地球的一部分的横截面图示出了不同类型的地层102、103、104，它们将包括本发明中的地震勘探数据。特别是，本领域的普通技术人员很快就会意识，本示例示出了共中点型道集，其中地震数据轨迹按表面几何形状分类，以模拟地球中的单个反射点。在这个例子中，来自几个炮点和接收器的数据可以组合成一个图像道集，或者根据要执行的分析的类型而单独地使用。尽管本示例可以说明平面反射器和相应的图像道集类别，但是也可以使用本领域中已知的其他类型或类别的图像道集，其选择可以取决于各种地球条件或事件的存在。
此外，如图1所示，来自多个入射点或震源104的地震能量将从不同地层之间的界面反射。然后，这些反射将被多个地震数据记录接收传感器105和井103捕获，每个传感器布置在彼此的不同位置偏移110处。因为所有入射点104和所有地震数据记录传感器105都放置在不同的偏移110处，勘探地震数据或轨迹(在本领域中也称为道集)将以各种入射角108来记录。入射点104在地球中生成向下的传输射线105，其通过记录传感器105由向上传输的反射来捕获。在该示例中，井位103用附接到井筒109上的现有钻井示出，使用本领域已知的技术可获得沿着井筒109的多次测量。该井筒109用于获得震源信息(如小波)以及测井数据，包括P波速度、S波速度、密度，以及其他地震数据。图1未示出的其他传感器布置在勘探区域101内，以便还捕获解释者和本领域的普通技术人员执行各种地球物理分析所需的层位数据信息。在本示例中，将根据现场记录来对道集进行分类，以检查振幅、信噪比、动校正、频率成分、相位和其他地震属性相对于入射角108、偏移测量110、方位角以及对数据处理和成像来说重要的且本领域的普通技术人员已知的其他几何属性的依赖性。
尽管入射点或炮点104在图1中表示为常见的共中点炮点几何形状，且炮点线大多是水平的，然而本领域的普通技术人员很快就会意识到，所述模式可以很容易地以其他方式表示，例如垂直的、对角线的或这三者的组合，这又会产生不同的炮点几何形状。然而，由于操作条件的原因，如图1所示的用于小波203、输入参数模型204和地震数据205的均匀采集的接收传感器105的均匀覆盖通常无法在整个勘探区域内实现。
参照图2a和2b，它们示出了本发明的优选实施例的概览的流程图。该方法始于检索阶段202，包括从勘探区域获得和准备数据和信息。特别地，从勘探区域检索三种不同类型的输入，包括：1)一组小波203，其源自在入射点104处产生的能量；2)输入参数模型集204，其还包括深度区间中的速度模型以及深度区间中的各向异性参数模型，两者均使用层析成像反演或全波形反演生成。本领域的技术人员很快就会认识到，这些模型可以包括深度区间的epsilon模型、深度区间的delta模型、深度区间的倾角模型、深度区间的方位角模型等；以及3)由接收传感器105直接从勘探区域101获得的地震数据205，其包括2D和3D地下结构数据，例如反射、剪切波和折射。所述检索阶段202由非暂时性计算机可读设备306执行，然后在步骤206处存储到内存资源304中，以供图3所示的数字高性能计算系统稍后进行处理和分析。
非暂时性计算机可读设备306然后从数据203、204和205的内存资源304中接收消息钩子，以启动检索过程207。非暂时性计算机可读设备306在步骤208处计算用于对震源小波203进行积分计算的预编程算法，使用公式：
Figure BDA0003139878800000151
该预偏移处理步骤208对于保证不失真的偏移频谱来说是必不可少的。该步骤不是为了产生真实的幅度RTM角度道集，而是与本实施例的拉普拉斯应用有关。在完成算法计算后，非暂时性计算机可读设备306在步骤209处生成新的震源小波，随后在步骤210处将其发送到内存资源304以进行存储。在成功存储新的震源小波209后，内存资源304向非暂时性计算机可读设备306发出消息钩子过程，以便计算机系统307通过用户界面在显示器309上显示消息，向终端用户指示该计算机系统307准备好了接收用户定义的有效截止频带值。此时，在步骤211，计算机实现方法的用户(例如本领域的普通技术人员)使用计算机系统设备输入的组合(例如键盘310和鼠标311)来在输入数据203、204和205的有效频带的范围内设置有效截止频带值。更准确地说，所述值的范围从0Hz到250Hz。一旦设置了所述值，用户将通过显示在显示器309上的消息对话框进行确认，该消息对话框将向非暂时性计算机可读设备306发送消息，以在步骤212将所设置的用户定义的有效截止频带存储到内存资源304。非暂时性计算机可读设备306然后在步骤213处从内存资源304中检索所设置的用户定义的有效截止频带、新的震源子波、输入参数204和地震数据205，并在步骤214开始计算逆时偏移算法。在处理逆时偏移计算的步骤214中，非暂时性计算机可读设备306在步骤215处在勘探区域的3D空间域方向x,y,z上生成原始的堆叠图像。所述生成的图像然后在步骤216由非暂时性计算机可读设备306存储到内存资源304。在成功完成存储216后，内存资源304向非暂时性计算机可读设备306发送消息，以启动从内存资源304中检索217所存储的原始堆叠图像，用于在步骤218中计算波数域x,y,z中的3D拉普拉斯算法和横向波数域x,y中的2D低通波数滤波算法。在步骤218执行的3D拉普拉斯算法是在空间域中的编程语言中定义的，其中P(x,y,z)是在步骤215生成的堆叠图像，包括以下公式：
Figure BDA0003139878800000152
此外，在步骤218，算法还被编程为使用
Figure BDA0003139878800000153
来计算波数域中的快速傅立叶变换FT[ΔP(x,y,z)]，其中
Figure BDA0003139878800000154
是P(x,y,z)的傅立叶变换；kx,ky,kz分别为x,y,z方向上的波数，由以下公式表示：
Figure BDA0003139878800000161
另一方面，在步骤218处，空间域中的2D低通波数滤波算法(也称为k-滤波)仅在横向(x,y)方向上进行，这与业内在所有方向上进行的其他方法不同。非暂时性计算机可读设备306然后在步骤219利用波数域x,y,z中的3D拉普拉斯算法和横向波数域x,y中的2D低通波数滤波算法来生成勘探区域的新图像，并在步骤220将其发送到存储资源304存储。在成功完成存储新图像的步骤之后，存储资源304向非暂时性计算机可读设备306发出信号，以在步骤221开始新图像的检索过程。在已经检索到图像之后，非暂时性计算机可读设备306在步骤222使用针对时间上的采样率的抗混叠算子来将图像从深度域转换为时间域，并给出可通过给定采样率来实现的高频地震道结果。该转换步骤222采用深度z方向，并使用来自输入参数204的速度模型和下述表达式将其转换到时域中：
Figure BDA0003139878800000162
非暂时性计算机可读设备306然后在步骤223生成针对深度z方向的时域中的经转换的图像。非暂时性计算机可读设备306然后在步骤224将所生成的转换到时域的深度z方向的图像存储到内存资源304。在成功地完成存储针对深度z方向的时域中的经转换图像的步骤224后，内存资源304向非暂时性计算机可读设备306发送信号，以在步骤225开始针对深度z方向的时域中的经转换图像的检索过程。非暂时性计算机可读设备306然后向计算机系统设备307发送信号，以在监视器309上显示编程的图形用户界面，其指示用户(通常是本领域的普通技术人员)确认已经设置了检索到的上述用户定义的有效截止频带的最大值(处于0Hz-250Hz之间范围的数字)，并且在步骤226继续将检索到的用户定义的有效截止频带的最大值设置为等于表达式fmax。非暂时性计算机可读设备306然后在步骤227使用步骤226的用户定义的有效截止频带的最大值fmax开始计算低通频率滤波。这样，计算机实现方法利用速度信息，将图像频率，准确地包含在设定的范围内并保持内容完整，从地震图像的浅层到其最深部分。因此，在步骤227期间，非暂时性计算机可读设备306自动实现频率内容对沿深度z方向的速度垂直变化的适应性，并从在步骤228处计算出的低通频率滤波中生成勘探区域的时域中的低通频率滤波图像。然后，非暂时性计算机可读设备306在步骤229处将所生成的勘探区域的时域中的低通频率滤波图像存储到内存资源304，并且在成功完成存储所生成的勘探区域的时域中的低通频率滤波图像的步骤229之后，内存资源304用信号通知非暂时性计算机可读设备306，以在步骤230开始针对所生成的勘探区域的时域中的低通频率滤波图像的检索过程。
在已经检索了图像的情况下，非暂时性计算机可读设备306在步骤231处使用在步骤222中使用的针对时间上的相同采样率的抗混叠算子来将图像从时间域转换为深度域，并给出在给定采样率下可达到的高频地震道结果。该转换步骤231采用深度z方向，并使用来自输入参数204的速度模型和下述表达式将其转换回深度域：
Figure BDA0003139878800000171
非暂时性计算机可读设备306然后在步骤232处生成勘探区域的3D空间域方向x、y和z中的最终堆叠图像。非暂时性计算机可读设备306然后在步骤232将勘探区域的3D空间域方向x、y和z中的最终堆叠图像存储到内存资源304。在成功地存储了勘探区域的3D空间域方向x、y和z中的最终堆叠图像之后，内存资源304用信号通知非暂时性计算机可读设备306，以在计算机系统307的显示监视器309上显示该计算机实现的逆时偏移方法，其使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，以生成增强的图像。
图3显示了数字高性能计算系统301的功能框图，该系统被编程为执行计算机实现的逆时偏移方法，其使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，从而生成增强的图像。该数字高性能计算系统301还包括(有线和/或无线的)内存资源304(其用于存储由接收传感器105使用无线传输系统305无线式传输的数据303)、非暂时性程序计算机可读存储器存储设备306，以及计算机系统设备307。
计算机系统设备307充当非暂时性程序计算机可读设备306的用户界面；输入、设置、选择和执行检索、计算、生成、调用、确定、转换和校正函数(消息钩子程序)的操作。所述计算机系统设备307连接到(有线和/或无线)非暂时性程序计算机可读设备306。计算机系统设备307还包括其他设备，如中央处理单元(CPU)308、显示器或监视器309、键盘310、鼠标311和打印机312。
数字高性能计算系统301具有固件、内核和软件，用于提供多个连接设备的连接和互操作性，例如用于存储数据的内存资源304、遥测系统305、非暂时性程序计算机可读存储设备306和计算机系统设备307。数字高性能计算系统301包括操作系统、一组消息挂钩过程和系统应用程序。
此外，由于性能始终是重要的问题，因此数字高性能计算系统设备301使用非暂时性程序计算机可读存储设备306以确保光束偏移步骤不会受到数字高性能计算的瓶颈system301 I/O，或任何网络通信。事实上，Apache Hadoop分布式文件系统和适当的数据压缩，以及根据数据的智能文件缓存将确保计算机实现的方法仅受内存/缓存速度和CPU计算能力的限制，而没有其他任何限制。
嵌入在数字高性能计算系统301中的操作系统可以是微软的“WINDOWS”操作系统、IBM公司的OS/2、UNIX、LINUX、Sun Microsystems或Apple操作系统，以及无数的嵌入式应用操作系统，例如可从Wind River,Inc.获得。
数字高性能计算系统301的消息钩子程序可以例如表示内存资源304、计算机系统设备307、非暂时性程序计算机设备306的操作或命令，其当前可能正在执行来自计算机实现的逆时偏移方法的特定步骤、过程或例程，该方法使用3D拉普拉斯算子、3D横向空间方向的低通波数滤波和1D时域中的低通频率滤波，用于去除低频和低波数噪声，从而生成增强的图像。
所述一组消息钩子程序可以首先由以下输入来启动：用户，例如输入用户定义的值或参数；计算机系统设备307的操作；在非暂时性程序计算机可读存储器存储设备306中的操作的处理；或者一旦某些数据已被内存资源304或非暂时性程序计算机可读存储器存储设备306所存储或检索到便自动执行。基于这些输入、过程或操作事件中的任何一个，内存资源304、非暂时性程序计算机可读存储器存储设备306或计算机系统设备307便生成数据包，该数据包被传送到系统计算机301，表示已发生的事件以及需要发生的事件。当系统计算机301接收到数据包时，它基于事件来将其转换为消息，并执行计算机实现方法的所需的步骤。该计算机实现方法包括一组消息钩子列表，其标识了一系列消息钩子程序。当操作系统收到消息后，将检查消息钩子列表，以确定是否有任何消息钩子程序已在操作系统中注册了自己。如果至少一个消息钩子程序已在操作系统中注册了自己，则操作系统会将消息传递给已注册的消息钩子程序，其会首先出现在列表中。被调用的消息钩子执行并向系统计算机301返回一个值，其指示系统计算机301将消息传递到下一个注册的消息钩子304、306或307。系统计算机301继续执行该操作，直到所有注册的消息钩子都通过为止，这通过识别磁干扰来表示该方法的完成。
根据本发明的优选实施例，对某些硬件和软件进行了详细描述，仅作为优选实施例，并不用于限制所公开实施例的实现结构。例如，虽然图1的接收系统设备的许多内部和外部组件。虽然已经描述了图3，但是本领域的普通技术人员将理解，这样的部件及其互连是众所周知的。此外，所公开的发明的某些方面可以体现在使用一个或多个接收系统、计算机系统设备或非暂时性计算机可读存储器设备执行的软件中。该技术的程序方面可以被认为是“产品”或“制品”，通常以可执行代码和/或相关数据的形式承载或体现在一种机器可读介质中。有形的非暂时性“存储”类型的介质和设备包括用于计算机、进程等的任何或所有存储器或其他存储器，或其相关模块，例如各种半导体存储器、磁带驱动器、磁盘驱动器、光盘或磁盘，以及之类的，可以随时为软件编程提供存储空间。
图4显示了该计算机实现的逆时偏移方法的图形用户界面，它显示在计算机系统307的监视器309上，示出了上述各种参数。
如本文所使用的术语“勘探区域”是指地质感兴趣的区域或体积，并且可以在任何测量尺度下与该区域或体积的几何形状、姿态和布置相关联。一个区域可能具有其中已发生的诸如折叠、断层、冷却、卸载和/或断裂的特征。
如本文所使用的术语“计算”包括各种各样的动作，包括计算、确定、处理、推导、勘探、查找(例如，在表、数据库或另一数据结构中查找)、确信等。它还可以包括接收(例如接收信息)、访问(例如访问内存中的数据)等。而且，“计算”可以包括解析、选择、选定、构建等。
如本文所用的术语“地表下”和“地下”是指在任何海拔高度上或在海拔范围内的任何一块陆地的顶面之下(无论是高于、低于还是处于海平面)，和/或低于任何地面的地表(无论是高于、低于还是处于海平面)。
除非另有明确说明，否则诸如“定义”、“创建”、“包括”、“代表”、“预分析”、“预定义”、“选择”、“构建”、“分配”、“创建”、“引入”、“消除”、“重新网格化”、“整合”、“发现”、“执行”、“预测”、“确定”、“输入”、“输出”、“识别”、“分析”、“使用”、“分发”、“干扰”、“增加”、“调整”、“合并”、“模拟”、“减少”、“分发”、“指定”、“提取”、“显示”、“执行”、“实现”和“管理”等术语可以指检索系统或其他电子设备的动作和过程，该系统或其他电子设备将某些电子设备的存储器(例如内存资源或非临时性计算机可读存储器)中的表示为物理量(电子，磁或光)的数据转换为存储器或传输设备或显示设备中的类似地表示为物理量的其他数据。
本文公开的实施例还涉及计算机实现的系统，用作用于执行本文所述的操作的检索系统的一部分。该系统可以被特别地构造用于所需目的，或者它可以包括由存储在内存资源或非暂时性计算机可读存储器中的计算机程序或代码选择性地激活或重新配置的通用计算机。这样，计算机程序或代码可以被存储或编码在计算机可读介质中，或在某种类型的传输介质上实现。计算机可读介质包括用于以诸如计算机之类的机器可读的形式存储或传输信息的任何介质或机构(“机器”和“计算机”在本文中可同义地使用)。作为一个非限制性示例，计算机可读介质可包括计算机可读存储介质(例如只读存储器ROM、随机存取存储器RAM、磁盘存储介质、光学存储介质、闪存设备等)。传输介质可以是双绞线、同轴电缆、光纤或其他一些合适的有线或无线传输介质，用于传输信号，例如电、光、声或其他形式的传播信号(例如载波、红外信号、数字信号等)。
本文所使用的接收系统或传感器105通常至少包括能够执行机器可读指令的硬件，以及用于执行产生期望结果的动作(通常是机器可读指令)的软件。另外，检索系统可以包括硬件和软件的混合体，以及计算机子系统。
硬件通常至少包括具有处理器功能的平台，例如客户端计算机(也称为服务器)和手持处理设备(例如智能电话、个人数字助理PDA或个人计算设备PCD)。此外，硬件可包括可以存储机器可读指令的任何物理设备，例如内存或其他数据存储器。其他形式的硬件包括硬件子系统，例如包括诸如调制解调器、调制解调器卡、端口和端口卡之类的传输设备。
软件包括存储在任何存储介质(例如RAM或ROM)中的任何机器代码，以及存储在其他设备(例如非暂时性计算机可读介质，例如外部硬盘驱动器或闪存)上的机器代码。软件可以包括源代码或目标代码，包含能够在客户端计算机、服务器计算机、远程桌面或终端中执行的任何指令集。
软件和硬件的组合也可以用于为所公开的发明的某些实施例提供增强的功能和性能。一个示例是直接将软件功能制造到硅芯片中。因此，应当理解，硬件和软件的组合也包括在检索系统的定义内，因此本发明将其设想为可能的等效结构和等效方法。
计算机可读介质或内存资源包括被动数据存储器，例如随机存取存储器(RAM)，以及半永久数据存储器，例如外部硬盘驱动器和外部数据库。另外，本发明的实施例可以体现在计算机的RAM中，以将标准计算机转换为新的特定计算机器。
数据结构是可以实现本发明的实施例的限定的数据组织。例如，数据结构可以提供数据的组织，或可执行代码的组织。数据信号可以跨非暂时性传输介质承载，并且跨各种数据结构存储和传输，因此可以用于传输本发明的实施例。
可以将系统计算机设计为可在任何特定的体系结构上工作。例如，该系统可以在高性能计算系统上执行，该高性能计算系统通常包括物理上连接的或通过局域网、客户端-服务器网络、广域网、互联网和其它便携式和无线设备及网络来连接的多个单台计算机的集合。
“输出设备”包括导致产生的直接行为，以及促进产生的任何间接行为。间接行为包括向用户提供软件，维护使用户能够通过其影响显示的网站，超链接至此类网站，或与执行此类直接或间接行为的实体合作。因此，用户可以单独操作或与第三方供应商合作，以便在显示设备上生成参考信号。显示设备可以作为输出设备包括在内，并且应适合于显示所需的信息，例如但不限于CRT监视器、LCD监视器、等离子设备、平板设备或打印机。显示设备可以包括已经通过使用旨在用于评估、校正和/或改善显示结果的任何常规软件进行校准的设备(例如，已经使用监视器校准软件进行了调整的彩色监视器)。作为在显示设备上显示参考图像的附加或替代，根据本发明的方法可以包括向对象提供参考图像。“提供参考图像”可以包括通过物理、电话或电子传递的方式将参考图像创建或分发给对象，通过网络提供对参考图像的访问，或者向被配置为在对象工作站上运行的对象或包含参考图像的计算机创建或分发软件。在一个示例中，提供参考图像可以涉及使对象能够经由打印机获得硬拷贝形式的参考图像。例如，信息、软件和/或指令可以被传输(例如，通过数据存储器或硬拷贝以电子或物理方式)和/或以其他方式可用(例如，通过网络)，以便于主体使用打印机来打印硬拷贝形式的参考图像。在这样的示例中，打印机可以是已经通过使用旨在用于评估、校正和/或改善打印结果的任何常规软件进行了校准的打印机(例如，已经使用颜色校正软件来调节的彩色打印机)。
一个数据库或多个数据库可以包含任何标准或专有数据库软件，例如Oracle，Microsoft Access，SyBase或DBase II。数据库可以具有字段、记录、数据和其他数据库元素，它们可以通过数据库专用软件进行关联。另外，可以映射数据。映射是将一个数据条目与另一个数据条目相关联的过程。例如，可以将包含在字符文件位置中的数据映射到第二个表中的字段。数据库的物理位置没有限制，并且数据库可以是分布式的。例如，数据库可能位于服务器的远程位置，并在单独的平台上运行。此外，可以在局域网和因特网的无线网络上访问数据库。
此外，模块、特征、属性、方法和其他方面可以被实现为软件、硬件、固件或其任何组合。在本发明的部件被实现为软件时，该部件都可以被实现为独立程序，作为较大程序的一部分，多个独立程序，静态或动态链接库，内核可加载模块，设备驱动程序，和/或计算机编程领域的技术人员现在或将来已知的所有其他方式。另外，本发明不限于在任何特定的操作系统或环境中的实现。
下面定义本文所使用的各种术语。对于在下面未对定义的在权利要求中使用的术语来说，应给予相关领域的技术人员所给出的反映在至少一个印刷出版物或授权专利中的最广泛可能的定义。
如本文所使用的置于第一实体和第二实体之间的“和/或”是指(1)第一实体、(2)第二实体以及(3)第一实体和第二实体中的一个。用“和/或”列出的多个元素应以相同的方式解释，即，如此连接的元素中的“一个或多个”。
另外，附图中的流程图和框图示出了根据本发明的各种实施例的系统、方法和计算机程序产品的可能实现的架构、功能和操作。还应注意，在一些替代实施方式中，方框中指出的功能可以不按图中指出的顺序发生。例如，取决于所涉及的功能，实际上可以基本上同时执行连续示出的两个方框，或者有时可以相反的顺序执行这些方框。还应注意，框图和/或流程图中的每个方框以及框图和/或流程图中的方框的组合可以由执行指定的硬件功能或动作的基于专用硬件的系统或者专用硬件和计算机指令的组合来实现。
尽管在前述说明书中已经针对本发明的某些优选实施例描述了本发明，并且出于说明的目的已经阐述了许多细节，但是本发明不应不当地限于已经出于说明性目的而阐述的前述内容。相反，在不背离如以下权利要求书所限定的本发明的真实范围的情况下，多种修改和替代实施例对于本领域技术人员将是显而易见的。另外应当理解，在本文的任何一个实施例中示出或描述的结构特征或方法步骤也可以在其他实施例中使用。

Claims (9)

1.一种计算机实现的逆时偏移方法，其中使用3D拉普拉斯算子、2D横向空间方向的低通波数滤波和1D时域的低通频率滤波，用于去除低频和低波数噪声，从而生成增强的图像，所述方法包括：
通过非暂时性计算机可读设备从具有空间域方向x、y和z的勘探区域中检索震源子波、输入参数模型集和地震数据；
将检索到的勘探区域的震源子波、输入参数模型集和地震数据存储到内存资源中；
通过所述非暂时性计算机可读设备从内存资源中检索所存储的震源子波、输入参数模型集和地震数据；
通过所述非暂时性计算机可读设备使用从内存资源中检索到的震源小波来计算积分算法；
通过所述非暂时性计算机可读设备从所计算的积分算法中生成新的震源小波；
通过所述非暂时性计算机可读设备将所生成的新的震源小波存储到内存资源中；
在所述非暂时性计算机可读设备上针对所检索到的地震数据来设置用户定义的有效截止频带；
通过所述非暂时性计算机可读设备来将所设置的用户定义的有效截止频带存储到内存资源中；
通过所述非暂时性计算机可读设备从存储资源中检索所存储的用户定义的有效截止频带、所存储的新的震源子波、输入参数模型和地震数据；
通过所述非暂时性计算机可读设备使用所检索到的用户定义的有效截止频带、所检索到的新的震源子波、所检索到的输入参数模型和所检索到的地震数据来计算逆时偏移；
通过所述非暂时性计算机可读设备根据所计算的逆时偏移来生成勘探区域的3D空间域方向x、y和z上的原始堆叠图像；
通过所述非暂时性计算机可读设备将所生成的勘探区域的3D空间域方向x、y和z上的原始堆叠图像存储到内存资源中；
通过所述非暂时性计算机可读设备从内存资源中检索所存储的原始堆叠图像；
通过所述非暂时性计算机可读设备使用所检索到的原始堆叠图像来计算波数域x,y,z中的3D拉普拉斯算法和横向波数域x,y中的2D低通波数滤波算法；
通过所述非暂时性计算机可读设备从所计算出的波数域x,y,z中的3D拉普拉斯算法和横向波数域x,y中的2D低通波数滤波算法中来生成勘探区域的新图像；
通过所述非暂时性计算机可读设备来将所生成的勘探区域的新图像存储到内存资源；
通过所述非暂时性计算机可读设备从内存资源中检索所存储的勘探区域的新图像；
通过所述非暂时性计算机可读设备使用所检索到的输入参数模型和时域中的抗混叠滤波器来将所检索到的勘探区域的新图像从深度域转换为时域；
通过所述非暂时性计算机可读设备生成勘探区域的时域中的经深度到时间转换的图像；
通过所述非暂时性计算机可读设备将所生成的勘探区域的时域中的经深度到时间转换的图像存储到内存资源中；
通过所述非暂时性计算机可读设备从内存资源中检索所存储的勘探区域的时域中的经深度到时间转换的图像；
通过所述非暂时性计算机可读设备将所检索到的用户定义的有效截止频带的最大值设置成等于表达式fmax；
通过所述非暂时性计算机可读设备使用所设置的用户定义的有效截止频带fmax的最大值来计算勘探区域的时域中的所检索到的经深度到时间转换的图像的低通频率滤波；
通过所述非暂时性计算机可读设备根据所计算的低通频率滤波来生成勘探区域的时域中的低通频率滤波图像；
通过所述非暂时性计算机可读设备将所生成的勘探区域的时域中的低通频率滤波图像存储到内存资源中；
通过所述非暂时性计算机可读设备从内存资源中检索所存储的勘探区域的时域中的低通频率滤波图像；
通过所述非暂时性计算机可读设备使用所检索到的输入参数模型和勘探区域的空间域方向z上的抗混叠滤波器来将所检索到的勘探区域的时域中的低通频率滤波图像从时域转换为深度域；
通过所述非暂时性计算机可读设备来生成勘探区域的3D空间域方向x、y和z上的最终堆叠图像；和
通过所述非暂时性计算机可读设备来将所生成的勘探区域的3D空间域方向x、y和z上的最终堆叠图像存储到内存资源。
2.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，通过非暂时性计算机可读设备从具有空间域方向x、y和z的勘探区域中检索震源子波、输入参数模型集和地震数据的步骤中的所述输入参数模型集还包括使用层析成像反演或全波形反演生成的深度区间的速度模型以及深度区间的各向异性参数模型。
3.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，所述非暂时性计算机可读设备还连接到计算机系统设备，用于将存储的步骤输出到显示器或打印设备。
4.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，通过所述非暂时性计算机可读设备使用从所述内存资源中检索到的震源小波来计算一组可执行指令的步骤还包括以下表达式：
Figure FDA0003139878790000031
5.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，在所述非暂时性计算机可读设备上针对所检索到的地震数据来设置用户定义的有效截止频带的步骤还包括处于0Hz到250Hz的值之间的用户定义的有效截止频带。
6.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，通过所述非暂时性计算机可读设备根据所计算的逆时偏移来生成勘探区域的3D空间域方向x、y和z上的原始堆叠图像的步骤进一步包括表达式：P(x,y,z)。
7.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，通过所述非暂时性计算机可读设备使用所检索到的原始堆叠图像来计算波数域x,y,z中的3D拉普拉斯算法和横向波数域x,y中的2D低通波数滤波算法的步骤还包括以下表达式：
Figure FDA0003139878790000032
8.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，使用所检索到的输入参数模型和时域中的抗混叠滤波器来将所检索到的勘探区域的新图像从深度域转换为时域的步骤进一步包括以下表达式：
Figure FDA0003139878790000041
9.根据权利要求1所述的计算机实现的逆时偏移方法，其特征在于，使用所检索到的输入参数模型和勘探区域的空间域方向z上的抗混叠滤波器来将所检索到的勘探区域的时域中的低通频率滤波图像从时域转换为深度域进一步包括表达式：
Figure FDA0003139878790000042
CN202110734440.XA 2020-06-30 2021-06-30 用于去除低频和低波数噪声以生成增强图像的方法和系统 Active CN113945982B (zh)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US16/917,571 2020-06-30
US16/917,571 US11333782B2 (en) 2020-06-30 2020-06-30 Computer-implemented method and system for removing low frequency and low wavenumber noises to generate an enhanced image

Publications (2)

Publication Number Publication Date
CN113945982A true CN113945982A (zh) 2022-01-18
CN113945982B CN113945982B (zh) 2023-07-25

Family

ID=79031800

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110734440.XA Active CN113945982B (zh) 2020-06-30 2021-06-30 用于去除低频和低波数噪声以生成增强图像的方法和系统

Country Status (2)

Country Link
US (1) US11333782B2 (zh)
CN (1) CN113945982B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20220283329A1 (en) * 2021-03-05 2022-09-08 Aramco Overseas Company B.V. Method and system for faster seismic imaging using machine learning
US11867856B2 (en) * 2021-06-03 2024-01-09 Saudi Arabian Oil Company Method and system for reflection-based travel time inversion using segment dynamic image warping
CN115755178B (zh) * 2023-01-06 2023-04-25 青岛欧谱赛斯海洋科技有限公司 一种基于积分地震子波的时间域全波形反演方法
CN116299696A (zh) * 2023-02-08 2023-06-23 中海石油(中国)有限公司深圳分公司 一种烃源岩泥质含量、toc、孔隙度、含水饱和度同时定量预测方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102269820A (zh) * 2010-06-01 2011-12-07 潜能恒信能源技术股份有限公司 一种基于gpu小存储量交错网格三维地震叠前逆时偏移成像方法
CN102608659A (zh) * 2012-03-27 2012-07-25 中国科学院地质与地球物理研究所 一种耦合透射系数的地震偏移方法
US20140133275A1 (en) * 2012-11-13 2014-05-15 Total E&P Usa, Inc. Process for Creating Image Gathers
CN103926623A (zh) * 2014-05-06 2014-07-16 王维红 一种压制逆时偏移低频噪音的方法
CN105510973A (zh) * 2014-09-23 2016-04-20 中国石油化工股份有限公司 一种用于压制逆时偏移成像噪音的拉普拉斯滤波方法
EP3076205A1 (en) * 2015-03-31 2016-10-05 CGG Services SA Method for survey data processing compensating for visco-acoustic effects in tilted transverse isotropy reverse time migration
US20160341837A1 (en) * 2015-05-20 2016-11-24 Charlie Jing Method of Removing Noise In Seismic Reverse-Time Migration
CN109975873A (zh) * 2019-04-23 2019-07-05 中国石油大学(华东) 一种逆时偏移成像去除低频噪音的方法及系统
MX2018010095A (es) * 2018-03-30 2019-10-01 Cgg Services Sas Metodos que emplean una inversion de forma de onda completa de tiempo de desplazamiento para produccion de imagenes de formaciones subsuperficiales con cuerpos de sal.

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8296069B2 (en) * 2008-10-06 2012-10-23 Bp Corporation North America Inc. Pseudo-analytical method for the solution of wave equations
US8619498B2 (en) * 2010-09-24 2013-12-31 CGGVeritas Services (U.S.) Inc. Device and method for calculating 3D angle gathers from reverse time migration
EP2771721A4 (en) 2011-10-28 2017-03-29 The Regents of The University of California Multiple-core computer processor for reverse time migration
WO2015104059A1 (en) * 2014-01-10 2015-07-16 Statoil Petroleum As Determining a component of a wave field
US20150276956A1 (en) * 2014-03-28 2015-10-01 Cgg Services Sa Wave-equation based processing and analysis of imaged seismic data
US9829592B2 (en) 2014-12-16 2017-11-28 Pgs Geophysical As Seismic imaging with visco-acoustic reverse-time migration using pseudo-analytical method

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102269820A (zh) * 2010-06-01 2011-12-07 潜能恒信能源技术股份有限公司 一种基于gpu小存储量交错网格三维地震叠前逆时偏移成像方法
CN102608659A (zh) * 2012-03-27 2012-07-25 中国科学院地质与地球物理研究所 一种耦合透射系数的地震偏移方法
US20140133275A1 (en) * 2012-11-13 2014-05-15 Total E&P Usa, Inc. Process for Creating Image Gathers
CN103926623A (zh) * 2014-05-06 2014-07-16 王维红 一种压制逆时偏移低频噪音的方法
CN105510973A (zh) * 2014-09-23 2016-04-20 中国石油化工股份有限公司 一种用于压制逆时偏移成像噪音的拉普拉斯滤波方法
EP3076205A1 (en) * 2015-03-31 2016-10-05 CGG Services SA Method for survey data processing compensating for visco-acoustic effects in tilted transverse isotropy reverse time migration
US20160341837A1 (en) * 2015-05-20 2016-11-24 Charlie Jing Method of Removing Noise In Seismic Reverse-Time Migration
MX2018010095A (es) * 2018-03-30 2019-10-01 Cgg Services Sas Metodos que emplean una inversion de forma de onda completa de tiempo de desplazamiento para produccion de imagenes de formaciones subsuperficiales con cuerpos de sal.
CN109975873A (zh) * 2019-04-23 2019-07-05 中国石油大学(华东) 一种逆时偏移成像去除低频噪音的方法及系统

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
柯璇等: "地震资料逆时偏移中的图形处理器加速算法", 《计算机系统应用》 *
柯璇等: "地震资料逆时偏移中的图形处理器加速算法", 《计算机系统应用》, vol. 22, no. 11, 30 November 2013 (2013-11-30), pages 115 - 118 *
牟海波等: "地震波逆时偏移中的层位校正与去噪方法", 《中国煤炭地质》 *
牟海波等: "地震波逆时偏移中的层位校正与去噪方法", 《中国煤炭地质》, vol. 28, no. 04, 30 April 2016 (2016-04-30), pages 71 - 75 *
郭念民等: "高阶拉普拉斯算子逆时偏移低频噪声去除方法", 《石油物探》 *
郭念民等: "高阶拉普拉斯算子逆时偏移低频噪声去除方法", 《石油物探》, vol. 52, no. 06, 30 November 2013 (2013-11-30), pages 642 - 649 *
陈可洋: "基于拉普拉斯算子的叠前逆时噪声压制方法", 《岩性油气藏》 *
陈可洋: "基于拉普拉斯算子的叠前逆时噪声压制方法", 《岩性油气藏》, vol. 23, no. 05, 31 October 2011 (2011-10-31), pages 88 - 95 *

Also Published As

Publication number Publication date
US20210405236A1 (en) 2021-12-30
CN113945982B (zh) 2023-07-25
US11333782B2 (en) 2022-05-17

Similar Documents

Publication Publication Date Title
Ikelle et al. Introduction to petroleum seismology
Prieux et al. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 1: Imaging compressional wave speed, density and attenuation
Virieux et al. An overview of full-waveform inversion in exploration geophysics
US9829592B2 (en) Seismic imaging with visco-acoustic reverse-time migration using pseudo-analytical method
US8352190B2 (en) Method for analyzing multiple geophysical data sets
CN113945982B (zh) 用于去除低频和低波数噪声以生成增强图像的方法和系统
US11493658B2 (en) Computer-implemented method and system employing nonlinear direct prestack seismic inversion for poisson impedance
US11143776B2 (en) Computer-implemented method and system for small cave recognition using seismic reflection data
KR20110057124A (ko) 지진 표면파들의 파형들을 사용하는 토양 특성들의 추정
WO2017035104A1 (en) Velocity model seismic static correction
Huang et al. Geological structure-guided initial model building for prestack AVO/AVA inversion
US10310117B2 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
CA2940406A1 (en) Characterizing a physical structure using a multidimensional noise model to attenuate noise data
CN113805237B (zh) 偏移陆地交叉排列地震的使用压缩感测模型的方法和系统
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
CN115877449B (zh) 用于在勘测区域内获得地下叠加图像的计算机实现方法
US11561312B2 (en) Mapping near-surface heterogeneities in a subterranean formation
Gras et al. Full-waveform inversion of short-offset, band-limited seismic data in the Alboran Basin (SE Iberia)
Rusmanugroho et al. 3D, 9C seismic modeling and inversion of Weyburn Field data
Zand et al. Integrated algorithm for high‐resolution crustal‐scale imaging using complementary OBS and streamer data
Kurzmann et al. Real data applications of seismic full waveform inversion
CN114114412B (zh) 使用各向异性参数来生成时间偏移图像道集的方法和系统
CN114442173B (zh) 预测和消除波束域中的多次波的计算机程序产品和方法
Bader et al. Source footprint elimination in full-waveform inversion by model extension: Application to elastic guided waves recorded by distributed acoustic sensing in unconventional reservoir
Consolvo Full-Waveform Inversion with Scaled-Sobolev Preconditioning Applied to Vibroseis Field 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
GR01 Patent grant
GR01 Patent grant