专利名称：用于对叠前地震数据成像的方法
技术领域：
本发明通常涉及地球物理学勘探领域。更具体的，本发明涉及地震数据处理领域。尤其是，本发明是用于对地震迁移的叠前地震数据进行成象的方法。
背景技术：
地震迁移是构造反射体表面的过程，该反射体表面根据记录的地震数据定义感兴趣的地下地层。因此，迁移的过程提供了在深度或时间中的地层的图像。试图说明与从地震源到地震接收器的地震能量的传输和反射相关的定位(运动学)和振幅(动力学)效应。尽管速度的垂直变化是最普遍的，横向变化也是会遇到的。这些速度变化是由一些原因引起的，其包括地层的分异致密，盐层侧面之上的升起，以及近端(页岩)和远端(含碳的沙)大陆架位置之间的沉积动态的变化。
许多年来已经认识到在地球物理处理工业中，可以通过叠前方式并且在深度域中实施地震迁移，而不是在堆栈之后，以便在具有复杂结构的区域中获得最佳的、堆栈图像。此外，叠前深度迁移提供了最佳制备用于后来AVO(相对偏移的振幅)分析的数据的优点。通过应用基尔霍夫(Kirchhoff)方法已经可以传统的实施叠前深度迁移。然而，由于计算硬件的新近发展和在有效外推器设计中的改进，开始使用以一维波动方程的解决方案为基础的方法。
已经在文献中很好证实了以波动方程为基础的方法提供了内含多路径的运动学优点，因此能够更精确的描绘复杂覆盖层下面的结构。然而，已经相当少的讨论可以通过基于波动方程的方法实现的动力学的优点。这并不意外，因为基尔霍夫叠前深度迁移相对于其以波动方程为基础的对手已经经历更长的发展。这种发展的一部分已经是各种因数或策略的发展，以在探测中解决几何偏差和空间不规则性。
叠后迁移常常不适用于地质复杂的区域，因为该迁移是在一稍后阶段，即在地震数据已经被NMO(正常迁出)或DMO(倾斜迁出)校正和CDP(共深度点)堆积影响之后执行的。一个叠前选择是普通拍摄收集或拍摄记录的深度迁移。该方法因此被称作是“拍摄记录”或“拍摄外形”迁移。拍摄记录迁移可以给出更精确的成像，更好的保持倾角，并且提供精确的叠前振幅信息。这些特性已经使拍摄记录迁移成为波动理论迁移的一个更普遍的方法。
在Reshef，Moshe和Kosloff，Dan著的“Migration of common shotgather”(Geophysics，51卷，第2期(February，1986)，第323-331页) 中描述了拍摄记录迁移。Reshef和Kosloff(1986)描述了用于迁移普通拍摄记录的三种方法。相比较于完全叠前迁移，拍摄记录迁移的精度和效率的讨论包含于两篇文献Jean-Paul Jeannot于1988年编的“Full prestack versus shot record mogrationpracticalaspect”，SEG(经济地质学家学会)Expanded Abstracts，第966-968页和C.P.A.Wapenaar and A.J.Berhout于1987年编著的“Fullprestack versus shot record migration”，SEG Expanded Abstract。
地震迁移包括两个步骤(1)波外推以及(2)成像。本申请解决在地震迁移中成像的第二步骤。在共同未批准的美国专利申请，“Method for Downward Extrapolation of Prestack Seismic Data”中描述了在地震迁移中波外推的相应方法，其和本专利申请具有同一发明者并指派于同一受让人。
通过测定来自单独的、平坦界面的用于反射响应的形式表达法可以获得用于精确拍摄记录迁移所必需的关键因素。然后可将拍摄记录迁移看成是用于平面波反射率反转所需的操作设置。在随后推广到多层情况之后，这些因素包括(1)源的定义，(2)使用保存垂直能量通量的振幅补偿运算在深度上外推波场相位，以及(3)实施在一些方面基于反射率定义的成像原理。
让上部半空间包含一个点，即具有坐标xs≡(xs，ys，zs)和频谱S(ω)的压缩地震源，并让它也包含地震接收器，比如具有坐标xr≡(xr，yr，zr)的水听器。正的向下选择深度，并且为了简化，假设zr＝zs。最后让界面的深度是z′。然后通过下面的式子给出由U标明的对于上行的(反射的)波场的压强的远场表达
U(xr,yr,zs,t)≈(i/c02)∫-∞∞S(ω)·exp(-iωt)∫-∞∞exp[ikk(xr-xx)]]]>·∫-∞∞exp[iky(yr-ys)]R[exp(-2iωξ0|z′-zs|)ωξ0]dkydkxdω---(1)]]>其中，c0和ξ0≡ξ(c0)是在上部半空间中的速度和垂直慢度，ω是辐射频率，kx和ky是水平波数，以及R是在反射器处的平面波反射率。在声学介质中，反射率R具有简单的形式R(kxω,kyω,z′)=(ρ1ξ1-ρ0ξ0)/(ρ1ξ1+ρ0ξ0)---(2)]]>其中，ρ0是上部半空间中的密度，以及ρ1、c1和ξ1≡ξ(c1)分别是在下部半空间中的密度、速率和垂直慢度。类似的，在反射物深度z′处通过D表示的下行波场具有形式D(kx,ky,z,ω)≈(i/c02)S(ω)exp[-i(kxxs+kyys)]exp(-iωξ0z)ωξ0---(3)]]>方程(3)提供了便利于理论分析的表示，但它对于用数字表示的应用不是最实用的形式。由于通常在(ω，x)域中执行拍摄记录迁移的外推，因此其更适宜代表在同一区域中的源。远离源，通过方程(3)的反变换给出波场D(x′,ω)=S(ω)exp(iω|x′-xs|/c0)c02|x′-xs|---(4)]]>其中，x′＝(x′，y′，z′)表示用于地下波场的观测点。
可以以至少两种方法使用方程(4)。第一，可在均匀覆盖层底部处指定波场，然后通过下面的、不同种类的结构向下连续。可替代的，在陆地设置中，利用局部速度在某一浅深度处指定波�。毓捶聪蛲馔浦帘砻妫缓笸ü旅娴摹⒉煌掷嗟慕橹氏蛳铝�。
这样，可以识别为平面波反射率R转换记录的地震数据所需要的步骤。首先，可能通过方程(4)，建立至反射物深度的源波场的模型。其次，上行的、记录的波场(U)外推至反射物深度z′。第三，通过下式计算反射率RR(kx,ky,z′,ω)=U(kx,ky,z′,ω)D(kz,ky,z′,ω)---(5)]]>并且然后回过来对空间域进行反变换。
在文献中已经报告了用于拍摄记录迁移的几个成像原理。观察五个通常的方法。运动学成像原理仅需要上行的、接收器波场U，而动力学成像原理依赖于在上行和下行、源波场U和D之间的某一函数关系。
通过下式给出普通的、运动学的方法R(x)=1nΣj=1nU(x,ωj)exp[iωjτ(x)],---(6)]]>其中，x＝(x，y，z)对应于成像点，而τ(x)是从射线轨迹确定的从源至成像点的单程时间。我们认为它能产生部分受射线理论的限制的结果。为了说明传播效应，必须从射线轨迹计算因数并以与基尔霍夫迁移相同的方式应用。对于该方法的进一步描述，例如，参见两篇文章Reshef，Moshe和Hosloff，Dan著的“Migration of common shotgather”，(Geophysics，51卷，第2期(1986年1月)，第324-331页)，orReshef，Moshe著的“Prestack depth imaging of three-dimensionalshot gather“，(Geophysics，56卷，第8期(1991年8月)第1158-1163页)。
通过下式给出在(x，ω)域中的平均相关性R(x)=1nΣj=1n[U(x,ωj)D*(x,ωj)],---(7)]]>其中，上标“*”代表复共轭。该方法将得到低解析(低频)的，但是加强的图像。然而，它没有得到反射率的有效测量。偏差效应的精确恢复需将反射率计算为源和接收器波场之间的比率的某种形式。
通过去卷积给出在(x，ω)域中的反射率的传统定义
R(x)1nΣj=1nU(x,ωj)D(x,ωj).---(8)]]>该成像条件非常接近于理论上的正确方法，其要求在波数中而不是在空间域中计算。该方法能够得到“浮点误差”，其中源照度非常弱。
通过下式给出在(x，ω)域中的标准化去卷积方法R(x)=1nΣj=1nU(x,ωj)D*(x,wj)|D(x,ωj)|,---(9)]]>其中，|D|＝(DD*)1/2。该运算法则相比于在方程(8)中的去卷积方法对噪声较不敏感。
在引入预白化之后，方程(9)变成R(x)=1nΣj=1nU(x,ωj)D*(x,wj)|D(x,ωj)|2+ϵ.---(10)]]>方程(10)是在工业中可能最通常使用的成像方法。该预白化项ε在存在噪声时稳定解。然而，它具有抑制几何偏差的精确恢复的不理想的副面效应。ε的最佳值高度依赖于数据并随每个拍摄记录变化并也随空间位置而变化。因此需要相当多的测试来发现单独的、“最佳的折衷值”。将尽可能小的对其进行选择，同时仍得到稳定、清楚看到的结果。对于参考方程(7)-(10)讨论的方法的概述，例如见1982年Jacobs，Allan著的“The prestack migration of profiles”，Ph.D.thesis，Stanford Univ.，SEP-34，和the classicbook，Claerbout，Jon，著的“Imaging the Earth`s Inventor”，Blackwell Scientific Publication，Ltd.，Oxford，UnitedKingdom，1985，1999。
Valenciano Alejandro，Biondi，Biondo，and Guitton，Antoine，于2002年著的，“Multidimensional imaging condition for shotprofile migration”，Stanford Exploration Project，ReportSEP-111，2002年6月10日，第71-81页阐明了一种用于拍摄记录迁移中的成象条件的最小平方反演。通过在位置x处计算反演时考虑空间临近位置处的数据，可以表示更稳定的(平滑的)最小平方反演。然而，在时域中而不是在频域中它们的反演全部被执行。在频域中，卷积和去卷积分别相应于乘法和除法。在时域中，去卷积包括带有时移操作数的矩阵构造。Valenciano等人(2002)形成的最小平方反演就是基于这些矩阵。
Valenciano等人(2002)也提出了用于上面的方程(10)中的预白化因数的简单变化。当结果UD*超过某一阀值时，可以加入二进位(0.或1.)加权以消除ε，而当UD*是不可接受的小时将它保持。因此，在清楚的成像区域中的结果将不会被预白化偏移，以及将在其它位置避免“浮点误差”。
因此，需要用于在拍摄记录迁移中对叠前地震数据进行成象的方法，其可以是完全自动的，并且恢复在不同种类的介质中的振幅。
发明内容
本发明是用于对叠前地震数据成象的方法。为地震数据中的每个频率计算单独反射率。然后通过单独反射率计算平均反射率。计算相对于频率的反射率组的方差。使用平均反射率为上行波场计算第二方差。接着使用反射率的方差和上行波场的方差计算空间变化的预白化因数。最后使用空间变化的预白化因数计算在每个位置处的单个反射率。
通过参考下面的详细描述和相关的附图，可以更容易的理解本发明和它的优点，其中图1是描述了用于对叠前地震数据成象的本发明的方法的实施例的处理步骤的流程；以及图2是描述了用于对叠前地震数据成象的本发明的方法的实施例的处理步骤的流程。
虽然结合优选实施例描述本发明，但可以理解本发明不受这些的限制。相反，本发明试图覆盖全部的替换、修改和等价物，只要其包括在通过附加权利要求限定的本发明的范围之内。
具体实施例方式
本发明是用于对叠前地震数据进行成象的方法。通过两个实施例描述本发明的方法。第一个实施例是迭代、非线性的反演以用于解决最佳、适应、空间变化的预白化因数。第二个实施例是最小平方的反演以确定在每个深度处为所有频率组确定最佳反射率。
在每个深度和每个横坐标处，施加动态(基于振幅的)成像原理以计算用于反射系数的估计。可以使用两种成像方法中的一个。第一成像方法通过在分母中加入预白化因数来计算所有“在下行之上的上行(upgoing over downgoing)”计算结果集的平均值，所述预白化因数在每个位置处自动适应噪声环境。该自适应是通过反复计算用于“数据”(接收器波场)和模型(反射率)的方差估计来实现的。第二成像方法在每个频率处结合使用“在下行之上的上行”计算的最小平方分析，从而为整组频率解决一个单独的反射率。
所述分析受限制于与一维波传播相关的通常假定，也就是忽略多次反射和传输损失。也假设材料特性在空间上平滑变化。在这些假设的限定范围中，拍摄记录迁移可以充分的恢复在不同种类的介质中的叠前振幅。
通过检测用于单独的界面的在(x，t)和(k，ω)域中的精确反射响应，可以获得用于拍摄记录迁移的成像。让x≡(x，y，z)是在具有速度c和与在具有时间相关性S(t)的xs处的点压缩源的无穷大的、均一的、声学的介质中的笛卡儿坐标。然后可将两维标量波方程写成▿2P(x,t)-(1/c2)∂2P(x,t)∂t2=-(4π/c2)δ(x-xs)S(t).---(11)]]>通过下式给出方程(11)的解P(x,t)=-S(t-|x-xs|/c)c2|x-xs|,---(12)]]>其中|x-xs|=(x-xs)2+(y-ys)2+(z-zs)2.]]>
通过下式给出下面的用于函数g(x，y，z，t)的正向缚立叶变换的协定g(kx,ky,z,ω)=∫-∞∞∫-∞∞∫0∞g(x,y,z,t)exp[iωt-kxx-kyy]dydxdt,]]>方程(12)中的解的Weyl积分形式是P(x,ω)=(i/2πc2)S(ω)∫-∞∞exp[ikx(x-xs)]]]>·∫-∞∞exp[iky(y-ys)]exp(-iωξ|z-zs|)ωξdkxdky,---(13)]]>其中ξ=1c2-kx2+ky2ω2.]]>对ω的方程(13)进行反向缚立叶变换将产生P(x,t)=(i/4π2c2)∫-∞∞S(ω)exp(-iωt)∫-∞∞exp[ikx(x-xs)]]]>·∫-∞∞exp[iky(y-ys)]exp(-iωξ|z-zs|)ωξdkxdkydω---(14)]]>方程(14)描述了由于在xs处的源，在空间位置x处观测到的直接获得的与时间相关的域。
下面，在深度z′处引入边界，并在(xr，yr，zr)处构造通过地震接收器，比如水听器记录的相关反射响应。为了方便，选择接收器的深度使其与用于源的深度相等，因此zr＝zs。基于“反射方法”，可以通过修改方程(14)的核心使其包括每个下行波成分的平面波反射系数来获得反射响应。对于在深度z′处的单独界面，反射率函数仅依赖于明显的慢度px＝kx/ω和py＝ky/ω，并且它不受频率的约束
R(kxω,kyω,z′)=(1ξ2-1ξ1)/(1ξ2+1ξ1),---(15)]]>其中ξ1和ξ2分别是在上部和下部半空间中的垂直慢度。因此通过下式给出在(xr，yr，zr＝zs)处的反射记录P(xr,yr,xs,t)=(i/4π2c2)∫-∞∞S(ω)exp(-iωt)∫-∞∞exp[ikx(xr-xs)]]]>·∫-∞∞exp[iky(yr-ys)]R(kxω,kyω,z′)exp(-2iωξ1{z′-zs|)ωξ1dkxdkydω,---(16)]]>这样，单独拍摄记录的迁移包括下列三个操作。第一，源波场从源的深度至成像深度向前外推，也就是模仿。第二，接受器波场从表面至成像深度反向外推。第三，施加成像原理，其可以是运动学的或是动力学的。使用这三个概念作为导向，可在单独的、一层响应中识别单独的波场。从等式(13)，通过下行波场D给出在反射器处源的前向外推场D(kx,ky,z′,ω)=(i/4π2c2)S(ω)exp[-i(kxxs+kyys)]exp(-iωξ1|z′-zs|)ωξ1.---(17)]]>因数“exp[-i(kxxs+kyys)]”来源于源的类似变数增值的、横向位置的傅里叶变换。由于从深度zs到深度z′的传播，因数“exp-(iωξ1|z′-zs|)”是相位的变化，并且因数“1/(ωξ1)”是用于每个平面波的振幅的垂直波数的权值。最后，因数“S(ω)”相应于源小波的傅里叶变换。
下面，使用方程(16)，在深度zr＝zs处记录的上行波场U被标示为U(kx,ky,zs,ω)=(i/4π2c2)S(ω)exp[-i(kxxs+kyys)]exp(-2iωξ1|z′-zs|)ωξ1R(ksω,kyω,z′).---(18)]]>为了根据反射物的深度编写U(Kx，Ky，Zs，W)，U必须使用从深度zs至z′的相反符号外推
U(kx,ky,z′,ω)=U(kx,ky,zs,ω)exp(iωξ1|z′-zs|)]]>=(i/4π2c2)S(ω)exp[-i(kxxs+kyys)]---(19)]]>exp(-iωξ1|z′-zs|)ωξ1R(kxω,kyω,z′).]]>比较方程(17)和(19)，在反射物处的下行和上行波场因而具有下面的关系U(kx,ky,z′,ω)=D(kx,ky,z′,ω)·R(kxω,kyω,z′),---(20)]]>其产生成像条件R(kzω,kyω,z′)=U(kx,ky,z′,ω)D(kx,ky,z′,ω).---(21)]]>如预期的，通过上行和下行波类型的比率这样给出平面波反射率。然而，该关系仅严格在明显的慢度或波数域中成立。通过kx和ky的反向傅里叶变换服从于π的因式范围之内，U(x，y，z′，ω)＝D(x，y，z′，ω)*R(x，y，z′，ω)， (22)其中“*”表示在坐标x和y之上的U和R之间的卷积。在方程(22)中已保持频率相关性，以便调整通过在许多频率处的测量提供的冗余。其通常实际上开始于渐近关系U(x，y，z′，ω)＝D(x，y，z′，ω)·R(x，y，z′，ω)， (23)其产生可替代的成像条件R(x,y,z′,ω)=U(x,y,z′,ω)D(x,y,z′,ω).---(24)]]>这样，对于同类介质，使用下面的运算法则可以完成振幅恢复拍摄记录迁移，其必须对每个深度z′进行重复。通过等式(17)，对于每个深度z′构造和外推源波场。然后，在通过kx和ky进行反向傅里叶变换以获得D(x，y，z′；ω)。傅里叶变换通过x和y对记录的波场U(x，y，z′；ω)进行傅立叶变换并反向外推至深度z’。然后，关于kx和ky进行反向傅里叶变换以获得D(x，y，z′；ω)。通过方程(24)计算R(x，y，z′)，其或者作为在全部频率之上的平均比率，或者作为使用预白化的关于标准化功率比率的平均数，或者如上面所述的通过最小平方。
方程(21)表示应该在波数域中实施成像。这是不便的，由于在空间区域中通�？梢曰竦猛ü邢薏罘滞馔频牟ǔ�。进一步的，方程(21)仅对平坦层有效，因为它内含零斜角的Snrll′法则的假设。
根据波阵面可以展开在空间区域中的反射响应。在高频处的第一次解成比例于下行波场的类似变数增值的行为。较高次的项成比例于海威得塞函数。该分析显示出对于一次项，在高频处，方程(24)是好的近似。我们仍然预料作为该近似的结果在横向结果中存在一些退化，其中在反射物特性中存在快速的横向变化，并且可能用于大的入射角度。
这样，渐近的，我们可以认为足以在(ω，x)域中计算比率。然后可以不依赖于倾角计算反射率。由于存在通过几个离散的频率提供的信息冗余，第一近似可以执行如下的平均R(x)=Σj=1nU(x,ωj)D*(x,ωj)|D(x,ωj)|2+ϵ(x)---(25)]]>其中，ε(x)是预白化项，其是趋于增加稳定性的小数。方程(25)相似于在上面的方程(10)中给出的反射率的功率标准化传统定义，但是具有可以空间变化的预白化项。
该空间变化预白化项ε(x)是必须的，因为对于一些频率，照度|D|可以依赖于噪声限制。在这种情况下，无限的反射率比率是非物理的并且可以异常大。在一些空间位置，可以保持这些限定用于大部分或甚至全部的频率成分，其需要ε随位置变化。
基于线性反演理论，可以从“模型”(δ2R)的方差的先前估计值和“数据”(δ2U)的测量方差评估空间变化的预白化项(ε(x))。不幸的是，较难猜测以前的评估。通过首先计算R的各个评估值以及用于每个频率的评估，我们可以得到用于δ2R的相当好的、初始评估，然后计算用于该分布的平均和标准偏差。然后可以与δ2U一起使用该产生的模型变化，以通过在方程(10)中设置ε(x)＝δ2U(x)/δ2R(x)获得最后评估。
这样，通过下面的自动程序可以处理在以经验为主地评估空间变化ε(x)过程中所涉及的困难。图1给出了表示用于对叠前地震数据进行成象的本发明的方法的第一实施例的处理步骤的流程。
首先，在步骤101，在叠前地震数据中选择多个n频率ωj。
在步骤102，选择初步的、恒定的预白化因数ε。初步的预白化因数不会在空间上变化，因此它在方程(25)中的应用，近似于在方程(10)中的较简单的反射率评估。
在步骤103，为在步骤101中选择出的每个频率ωj计算反射率R。优选的，利用在步骤102中选择出的初步的、恒定的预白化因数使用方程(25)计算反射率。
在步骤104，根据在步骤103中计算的反射率的分布计算平均反射率<R(x)>。利用如本领域熟知的任何合适的方法计算平均反射率。
在步骤105，根据在步骤103中计算的反射率的分布计算用于反射率的方差δ2R(x)。优选的，利用下面的方程计算用于反射率的方差。
σR2(x)=1n-1Σj=1n[U(x,ωj)D(x,ωj)-(R(x))]2.---(26)]]>在步骤106，利用在步骤104中计算的平均反射率<R(x)>计算上行波场的方差δ2U(x)。优选的，利用下面的方程计算上行波场的方差。
σU2(x)=1n-1Σj=1n[(x,ωj)-(R(x))·D(x,ωj)]2.---(27)]]>在步骤107，利用在步骤105中计算的反射率的方差和在步骤106中用于计算的上行波场的方差计算空间变化的预白化因数ε(x)。
优选的，通过下面的方程计算预白化因数ϵ(x)=σU2(x)σR2(x).---(28)]]>在步骤108，计算反射率R(x)。通过再次对方程(25)施加在步骤107中计算的空间变化的预白化因数计算反射率，因此，通过下面的方程优选计算反射率R(x)=1nΣj=1nU(x,ωj)D*(x,ωj)|D(x,ωj)|2+σU2(x)σR2(x).---(29)]]>在步骤109，确定在步骤10中计算的反射率R(x)是否是令人满意的。如果答案是否定的，则反射率不是令人满意的，然后处理返回至步骤105以计算另一个反射率值。如果答案是肯定的，则反射率是令人满意的，然后处理中止于步骤110。
一个以上的迭代是必需的，其依赖于数据的特征。理想的，通过循环方程(26)至(29)可以实施几次迭代。但是，为了速度，该优选实施例仅实施了一次迭代。
也可以构造本发明的方法的第二实施例。在方程(6)到(10)中描述的成像原理全部表示关于每个频率ωj的反射率R的单独估计的某一平均数。每个估计是具有零截距的直线的斜率，其涉及上行至下行波场。因此，平均反射率是斜率的平均值。作为替代，本发明的方法的第二实施例是用于单个斜率的最小平方解，其使受到线具有零截距限制的剩余误差的平方和最小化。现在，将D和U作为关于频率的向量对待。那么将用于单个参数R的反演问题采用下面的形式U(x，ωj)＝D(x，ωj)·R(x).(30)方程(30)具有最小平方解
R(x)=DH(x)·U(x)DH(x)·D(x),---(31)]]>其中，“H”代表转置矩阵(赫尔曼变化)的复共轭。增加预白化因数ε，在方程(31)中的反演变成R(x)=DH(x)·U(x)DH(x)·D(x)+ϵ.---(32)]]>使用方程(30)的向量表示，方程(32)可以再次写成R(x)=1nΣj=1nD*(x,ωj)U(x,ωj)1nΣj=1nD*(x,ωj)D(x,ωj)+ϵ.---(33)]]>在该公式中，不必使用空间变化的预白化因数ε(x)。可将方程(33)应用于空间或波数域。由于R的这个新的估计被计算作为和的比率，而不是比率的和，所以即使在缺少预白化因数ε的情况下，我们将遇到更少的调节问题。
这样，对于不同种类的介质，可以使用下面的方法完成振幅恢复拍摄记录迁移。图2给出了描述用于成像叠前地震数据的本发明的方法的实施例的处理步骤的流程。
首先在步骤201中，为叠前数据组选择多个深度z′。从本领域公知的任何合适的源选择深度，比如用于地震数据组的可用速度模型。
在步骤202中，从在步骤201中选择出的多个深度选择出深度z′。以系统化的方式选择深度，所述方式开始于上部深度，并通过地震数据组向下逐渐变深。
在步骤203中，构造源波场D(x，y，z′；ω)并外推至在步骤202′中选择的深度z′。可以例如，通过方程(17)获得最后的源波场D(x，y，z′；ω)。
在步骤204中，记录的波场U(x，y，z′；ω)关于x和y进行傅立叶变换并且反向外推至深度z′。可以例如，通过方程(18)和(19)获得最后的波场U(x，y，z′；ω)。
在步骤205中，在深度z′处计算反射率R(x，y，z′)。可以例如，通过方程(21)计算反射率。可以将反射率计算作为关于全部频率的比率的平均数，作为关于使用预白化的标准化功率比值的平均，或者通过最小平方计算反射率。优选地，通过上面描述的本发明的两个实施例中的一个计算反射率。本发明的第一实施例，参考上面的方程(25)-(29)讨论了带有空间变化预白化因数的标准化功率比值。参考上面的方程(33)讨论了本发明的第二实施例，即最小平方的方法。
在步骤206中，确定任何深度z′是否保持为在步骤201选择出的多个深度中选择出的。如果答案是肯定的，则保持深度z′，然后处理返回至步骤202以选择下一个深度。如果答案是否定的，则不保持深度z′，然后处理中止于步骤207。
可以理解的是前面的描述仅仅是本发明的具体实施例的详细说明，在不脱离本发明的范围得情况下，根据公开的内容可对披露的实施例产生大量的变化、修改和替代，只要其。因此，前述描述不是用于限制本发明的范围的。相反，仅通过附加权利要求和它们的等价物确定本发明的范围。
权利要求
1.用于对叠前地震数据成象的方法，包括为地震数据中的每个频率计算单独的反射率；根据所述单独的反射率计算平均反射率；为所述单独反射率计算方差；使用平均反射率为上行波场计算方差；使用反射率的方差和上行波场的方差，计算空间变化的预白化因数；以及使用空间变化的预白化因数计算反射率。
2.权利要求1的方法，其中为上行波场计算方差的步骤包括应用下面的方程σU2(x)=1n-1Σj=1n[U(x,ωj)-⟨R(x)⟩·D(x,ωj)]2,]]>其中，δ2U(x)是上行波场的方差，x是空间位置，n是频率ωj的个数，U(x，ωj)是上行波�。�<R(x)>是平均反射率，以及D(x，ωj)是下行波场。
3.权利要求1的方法，其中计算空间变化的预白化因数的步骤包括应用下面的方程ϵ(x)=σU2(x)σR2(x),]]>其中，ε(x)是空间变化的预白化因数，δ2U是上行波场的方差，以及δ2R(x)是反射率的方差。
4.权利要求1的方法，其中使用空间变化的预白化因数计算反射率的步骤包括应用下面的方程R(x)=1nΣj=1nU(x,ωj)D*(x,ωj)|D(x,ωj)|2+σU2(x)σR2(x),]]>其中，R(x)是反射率，x是空间位置，n是频率ωj的个数，U(x，ωj)是上行波场，D(x，ωj)是下行波�。�2U是上行波场的方差，以及δ2R(x)是反射率的方差。
5.一种用于对叠前地震数据成象的方法，其中最小平方法包括应用下面的等式R(x)=1nΣj=1nD*(x,ωj)U(x,ωj)1nΣj=1nD*(x,ωj)D(x,ωj)+ϵ,]]>其中，，R(x)是反射率，x是空间位置，n是频率ωj的个数，U(x，ωj)是上行波�。珼(x，ωj)是下行波�。约唉攀窃ぐ谆蚴�。
全文摘要
通过为地震数据中的每个频率计算单独反射率而对叠前地震数据进行成象。然后通过单独反射率计算平均反射率。计算相对于频率的反射率组的方差。使用平均反射率为上行波场计算第二方差。接着使用反射率的方差和上行波场的方差计算空间变化的预白化因数。使用空间变化的预白化因数计算在每个位置处的反射率。
文档编号G01V1/28GK1609634SQ20041008705
公开日2005年4月27日 申请日期2004年10月22日 优先权日2003年10月23日
发明者S·M·凯利 申请人:Pgs美洲公司
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
