专利名称：三维双向声波动方程叠前成像系统和方法
相关申请的交叉引用无适用。
关于联合发起的研究或开发的声明无适用。
背景技术：
在地震探测中，通过将声能传进地球的近地表面并且探测自地下地球构成的不同层之间的边界反射的声能来获得地震数据。在位于边界相对两侧上的层之间存在有声阻抗差异时声能被反射。表示已探测声能的信号被解释以此推断地下地球构造的结构和组成。
在海洋地震探测中，(在水体中实施的地震探测)地震能量来源(如气枪或气枪阵列)通常用来将声能传进地球。气枪或气枪阵列在水中选定的深度处被驱动，通常同时气枪或气枪阵列被地震勘探船所牵引。相同或不同地震勘探船在水中还拖着一条或多条被称为“拖缆”的地震传感器电缆。通常，拖缆在船后沿着其中拖缆被牵引的方向伸展。一般，拖缆包括多个压力传感器，通常为水听器，它们以一定距离的间隔被布置在电缆上，其沿着电缆的位置是已知的。水听器是生成与水压或水中压力的时间梯度(dp/dt)相对应的光或电信号的传感器。牵引一条或多条拖缆的船舶通常包括记录仪器，用来制成由水听器响应已探测的声能而生成的信号的、根据时间索引的记录。如前面所解释的，信号的记录被处理以此推断实施地震勘探的位置下面的地球构成的结构和组成。在这里，所描述的本发明是一种用于处理数据的改进方法，并且具体地说是一种用于地震数据的深度偏移的改进方法。
在考虑了下面的详细说明以及附图之后可以获得对各个不同的公开实施例的更好的理解，其中图1示出了说明性的海洋地震勘探环境；图2示出了说明性的地震勘探记录系统；图3示出了一组说明性的记录信号；图4示出了说明性的三维数据体；图5示出了说明性的成像系统；图6示出了用于说明性的成像方法的数据流；以及图7示出了用于说明性的成像方法的流程图。
虽然本发明允许各种修改和替换形式，但是其中的特定实施例通过举例被示出在附图中并且在这里将被详细描述。然而，应当理解，附图以及对此的详细描述不是用来将本发明限定在所公开的特定形式上，而是相反，本发明将覆盖属于如所附权利要求所界定的本发明范围内的所有的修改、等同物和替换。
具体实施例方式
下文描述了用于由叠前地震共炮集导出图像信息的系统。利用高阶空间有限差分将双向声波动方程用来产生自源位置随时间向前传播以及自接收器位置随时间向后传播的波场。成像标准用于源和接收器波场以此获得分图像，所述分图像在空间上被相加以此获得最终的图像。本发明的某些实施例使用声波动方程的最大时间步进和时间抽选，以此加快计算并减少存储需求以及有效的任务分配策略，从而在多个处理器当中并行计算。
地震探测依照常规在陆地和海上进行。在海上，地震勘探船在船后配备拖缆，如图1所示。当船舶向前移动时(在箭头102的方向上)，每条拖缆110拖在船100后，并且每条拖缆包括多个接收器114。每条拖缆110还可包括可编程转向器118和可编程深度控制器，其将拖缆抽出到达距船舶路径的操作偏距以及向下到达操作深度。
拖缆110可以长达几英里，并且通常被构造成若干部分，每部分的长度为25至100米，所述的若干部分包括若干组多达35个或以上的接收器。每条拖缆110包括电或光缆，用于互连接收器114和船舶100上的地震仪器。数据在接收器114附近被数字化并且通常以每秒7(或以上)兆比特的数据速率通过电缆线路被传输到船舶100。
如图1所示，地震勘探船100还牵引源112。源112可以是脉冲源或振动源。在海洋地震学中使用的接收器114被总称为水听器，并且通常利用压电换能器来构造。各种适当类型的水听器是可用的，如圆盘水听器和圆柱水听器。船上的仪器控制源和接收器的操作并记录所获取的数据。
地震勘探提供数据用于评估地下结构的位置，如结构106，其位于海底108下。通过估计至地下结构的各种距离可确定结构的几何形状或地形。记录的地震数据的某些幅度和其他属性表示油床或气储集层的存在。
为了确定到地下结构106的距离，源112发射地震波116，地震波116反射离开地下结构106(和其他地下结构)。反射的波被接收器114感测。通过确定地震波116从源112传播到地下结构106到接收器114的传播时间，可获取到地下结构106的距离的估计值。
图2示出了说明性的地震勘探记录系统，该地震勘探记录系统具有耦合至总线202的接收器114以此将数字信号传递至勘测船100上的数据记录电路206。位置信息和其他参数传感器204也耦合至数据记录电路206以使数据记录电路能够存储另外的信息，比如对解释纪录的数据有用的阵列取向信息。
通用数字数据处理系统208被示出耦合至数据记录电路206，并且还被示出通过总线202耦合至定位装置210和地震源112。处理系统208配置记录电路206、定位装置210和地震源112的操作。记录电路206在非易失性存储介质(如光盘或磁盘的存储阵列)上获得来自接收器114的高速数据流。定位装置210(包括可编程转向器和深度控制器)控制接收器114和/或源112的位置。
图2的地震纪录系统可包括在这里没有明确示出的其他部件。例如，处理系统208可包括具有图像显示的用户界面和键盘或接受用户输入的其他方法。每条拖缆110可具有独立总线202，用于耦合到数据记录电路。处理系统208还可包括网络接口，用于将存储的地震勘探数据传递至具有强大计算资源的中央计算设施用于处理地震勘探数据。
图3示出了被接收器114探测并取样的说明性的地震信号。信号表示作为时间的函数的地震波能量的某些测量(如位移、速度、加速度、压力)，并且它们在可编程的取样速率下以高分辨率(如24位)被数字化。这样的信号可以不同的方式来分组，并且在如此分组时，它们被称为“集”。例如，“炮集”是为地震源的单次开火纪录的信号组。
纪录的地震勘探数据在保持图3的格式时是没有多少用处的。尽管能够在图中并排绘制揭示大尺度地下结构的各种纪录波形，但是这样的结构是失真的并且甚至看不到更精细的结构。这种技术扩展到三维是难以使用的。图4说明了优选的数据格式。
图4示出了“数据立方体”，即数据值的三维阵列。(不管术语“立方体”的用法，数据立方体的维度不必相等。)数据立方体表示整个勘探区域的某些地震属性。三维阵列包含大小均匀的单元，每个单元具有表示该单元地震属性的数据值。可以表示各种地震属性，并且在某些实施例中，每个单元具有表示多种地震属性的多个数据值。适当的地震属性的实例包括反射率、声阻抗、声速和密度。数据立方体格式有助于计算分析和直观描绘，并且因为这个原因，数据立方体可以称为勘探区域的“三维图像”。
图5示出了处理系统500，用于将存储的接收器信号转换成勘探区域的三维图像。个人工作站502通过局域网(LAN)504耦合至一个或多个多处理器计算机506，其又通过存储区域网(SAN)508耦合至一个或多个共享存储器510。个人工作站502用作到达处理系统的用户接口，使得用户能够将勘探数据装载进系统，从系统中取回图像数据并且配置和监控处理系统的操作。个人工作站502可采用具有图形显示和键盘的台式计算机的形式，所述的图像显示以图形显示勘探区域的三维图像和勘探数据，所述键盘使得用户能够移动文件并执行处理软件。
LAN504提供了多处理器计算机506之间的以及同个人工作站502的高速通信。LAN504可采用以太网的形式。
多处理器计算机506提供了并行处理能力，以此适合使地震接收器信号能够快速转换成勘探区域图像。每个计算机506包括多个处理器512、分布式存储器514、内部总线516、SAN接口518和LAN接口520。每个处理器512处理输入数据的分配部分以此产生地震勘探区域的分图像。与每个处理器512相关联的是分布式存储模块514，分布式存储模块514存储处理器使用的转换软件和工作数据集。内部总线516提供处理器间通信以及通过相应接口518、520与SAN或LAN网络的通信。在不同计算机506上的处理器之间的通信可以通过LAN504来提供。
SAN508提供了对共享存储装置510的高速存取。SAN508可采用光纤通道或无限宽频网络的形式。共享存储器510可以是将磁盘介质用于非易失性数据存储的大的独立信息存储器。为了改进数据存取速度和可靠性，共享存储器510可以配置成冗余磁盘阵列。
图6示出了用于通过处理系统500实现的某些转换过程实施例的数据流。共享存储器510最初存储速度立方体602和炮集604。炮集604是若干组来自地震勘探的存储的接收信号。速度立体方602是声速的数据立方体，并且它可以利用已知技术从炮集604中导出。参见比如Jon F.Claerbout的Fundamentals of Geophysical DataProcessingp.246-56，其内容通过引用被结合在此。
每个处理器512从共享存储器510中取回炮集和速度立方体的其中之一。以源位置和测得的或模拟的源波形为起点，每个处理器512向前传播源波形以此生成作为时间的函数的源波场608。源波场608可以均匀间隔的时间间隔被保存在存储器510中。接着，以接收器位置和接收的波信号为起点，每个处理器512向后传播接收的波形以此生成作为时间的函数的接收器波场608。至多，每个处理器512可预计在滑动时间窗口中保存接收器波场。当接收器波场608被产生时，它们可以利用成像标准与相应的源波场组合以此生成分图像。每个处理器512组合分图像以此建立图像区域610，并且重叠的图像区域在空间上被相加以此建立勘探区域的完整图像，该图像被存进共享存储器510中。
图7示出了说明性的转换方法实施例的流程图。从块702开始，处理系统500获取地震勘探数据。勘探数据包含纪录的接收器波形，其可以针对单个接收器振动来校正，但是在其他的方面它们是原始的、未偏移的叠前数据。如果数据还未组织成炮集604，系统500可实施这种数据分组。
在块704，处理系统500获取以前已经在别处利用本领域普通技术人员已知的方法确定的速度数据立方体602。在块706，处理系统500确定空间栅格大小、时间步长和存储抽选因子。这些值的每一个可以作为输入参数由用户来提供，或者依照下列方程由处理系统来计算。空间栅格大小是对波场数据立方体中单元尺度Δx，Δy，Δz的限制。其被计算如下max(Δx，Δy，Δz)＝cmin/(αfmax)(1)其中cmin是图像空间中的最小声速，fmax是源小波的(或输入数据中的)最大频率，以及α是最大频率波的每波长的样本数。每fmax波长的样本数α是处理系统500被配置成可使用的波传播方程的逼近阶的函数。对于2nd阶空间差分逼近，α应当是10。对于4nd阶，α是4，并且对于更高阶空间差分逼近，3是α的充分值。(参见M.A.Dablain的“The application of high-order differencing to scalar waveequation，”GEOPHYSICS，Soc.Of Expl.Geophys.，v.51，p.54-66，1986)。
模拟时间步长Δt限制接着可确定为Δt≤μmax(Δx，Δy，Δz)/cmax(2)其中cmax是图像空间中的最大声速，并且μ是取决于空间差分滤波器的时间步长因子。时间步长因子可按下式来确定
μ=2Σp=-n/2n/2(|wpx|+|wpy|+|wpz|)]]>(3)其中n是空间差分滤波器的阶，并且wpx，wpy，wpz是x，y和z方向的空间滤波器系数。(空间滤波器系数将在下面作进一步讨论。)已经确定了模拟时间步长Δt，处理系统500接着确定存储抽选因子(或许取回它作为用户指定的参数)。为了减少存储需求并加快计算，处理系统未保存源波场的每个时间步长。相反，每kth时间步长(如每10th时间步长)被存到磁盘中。在某些备选实施例中，还可使用空间抽选或源波场压缩来进一步减少存储需求。
在块708，处理系统500确定处理器任务分配。在某些实施例中，每个炮集被分配给相应处理器，其中每个处理器可能处理多个炮集。在其他实施例中，给定炮集的波传播在多个处理器之间进行分配，其根据需要交换信息以此利用区域分解方法完成源和接收器波场。(参见A.Villarreal和J.A.Scales的“Distributed threedimensional finite-difference modeling of wave propagation in acousticmedia，”COMPUTERS IN PHYSICS，American Inst.Of Physics，v.11，p.388-99，1997。)块710描绘了多个处理器可并行执行的计算。每个处理器实施由块712-722表示的、针对其分配任务的每一个的操作。在块712，每个处理器512获得特定任务参数。例如，当针对每个炮集的处理被分配给单个处理器时，块712表示炮集数据的取回，包括源位置和接收器位置，由此可确定模拟区域的位置和尺度。如果每个处理器只处理与给定炮集相关联的波场部分，每个处理器另外确定哪一部分将要处理以及处理该炮集的其他处理器。
在块714，每个处理器512利用向前传播生成作为时间的函数的源波场。声波动方程由下式给出1c2∂2u∂t2=▿2u+s=∂2u∂x2+∂2u∂y2+∂2u∂z2+s---(4)]]>
其中u＝u(x，y，z，t)是压力场，c＝c(x，y，z)是由速度立方体给出的速度场，以及s＝s(x，y，z，t)表示来自地震源的激励信号。拉普拉斯算子2u可以在矩形栅格上逼近nth阶(n是偶整数)(▿2ux,y,zt)n=Σl=-n/2n/2wl(ux+l,y,ztΔx2+ux,y+l,ztΔy2+ux,y,z+ltΔz2)---(5)]]>其中ux，y，zt是坐标x，y，z和时间t下压力场的值，以及wl＝w-l是空间有限差分权重。空间有限差分权重利用级数展开法来确定。计算高阶空间有限差分权重的通用公式可以按照下列符号表达式给出{w-no/2,···,w-1,w0,w1,···,wno/2}=Coefficients[Series[sn/2(lnsΔx)2,{ons,arounds=l,ton}]]---(6)]]>其中n是空间差分的阶并且它是偶数，系数[T(s)]是表示多项式T(s)的系数集合的表达式，级数[A，{B，C，D}]是变量B在C点处的表达式A的D阶泰勒级数展开。在表达式A中，参数s被定义为eiωΔx。
利用依照方程(2)选择的时间步长，二阶有限时间差分对于方程(4)中的时间导数是适当的。因此，方程(4)的离散化型式变为ux,y,zt+1=2ux,y,zt-ux,y,zt-1+Δt2cx,y,z2((▿2ux,y,zt)n+sx,y,zt)---(7)]]>其中cx，y，z是坐标x，y，z处的声速值，并且sx，y，zt是坐标x，y，z和时间t下的地震源信号。用于向前传播的地震源信号sx，y，zt可以是来自地震源112的记录的或模拟的信号。另一方面，还可利用简单源小波。例如，可使用脉冲函数、限带零相位雷克小波或另一种类型的零相位小波。
方程(7)被用于向前的波传播。对于向后的波传播(块716)，离散化方程是ux,y,zt-1=2ux,y,zt-ux,y,zt-1+Δt2cx,y,z2(▿2ux,y,zt)n---(8)]]>在向后的传播中，接收的信号被用来设置适当的压力场值，与其中地震源增加适当的压力场值的向前的传播不同。
在实施(离散化的)波动方程的向前和向后传播时，处理器512强迫服从各种边界条件。为了处理表面相关的多个反射，自由表面边界条件在海洋表面位置处被设置。模拟空间的其余边界在真实世界中是不存在的，并且因此，处理器512使用“吸收边界”来抑制不想要的、自这些边缘的反射。吸收边界包括具有锥度的空间对模拟空间的填充，从而导致额外的有限差分层起到象吸收海绵的作用。声波场可沿余弦锥度或指数锥度而行，从模拟空间内的波场的100％处开始并且下降到大约60-80％。
返回到图7，在块714，每个处理器512实施离散化波动方程的向前传播。源波场是压力场的时间级数。压力场以规律间隔的时间步长被保存在本地磁盘或共享存储510上用于将来参考。
在块715，每个处理器512将用于向后传播的时间索引设为等于来自块714的最后的时间索引。在块716，每个处理器512实施离散化波动方程的向后传播以此生成接收器波场。以与存储的源波场相对应的时间步长，接收器波场可以被保存用于将来参考。另一方面，当生成接收器波场时，还可以通过将成像条件施加于接收器波场和相应的源波场来避免接收器波场的保存。
在块718，每个处理器52将成像条件施加于具有相应的时间步长的源和接收器波场。在某些实施例中，成像条件是在时间上积分时为每个共炮集生成下面的分图像的两个波场的零滞后互相关Px,y,z=Σmfx,y,zmkΔtbx,y,zmkΔtΣxΣyΣz|fx,y,zmkΔt|2+σ---(9)]]>其中kΔt是用于保存的源波场的时间步长(在通过k的抽选之后)，fx，y，zt和bx，y，zt是坐标x，y，z和时间t下的用于源和接收器波场的压力场值，以及σ是选择用来防止消散分母导致数值不稳定性的小的恒定值。因为方程(9)上的分母补偿了源波场的强度，它可被认为是“照明校正”。分图像揭示了向前和向后传播波场随时间处处匹配的反射器。在块719，进行检查以此了解时间索引是否已经达到零，并且如果未达到零，则每个处理器512返回到块716。
在块720，每个处理器512将分图像存进整个勘探区域内的适当位置，在空间上对处处存在重叠的具有以前存储值的分图像求和。在地震勘探期间，当船舶(和接收器阵列)向前移动时，地震源被重复点火。空间源点火间隔被设计成可提供连续炮集之间的大的重叠，从而使得整个勘探区域的每一部分能够根据多个分图像来确定。
在块722，每个处理器512确定另外的任务(如炮集)是否已经被分配给该处理器。如果另外的任务存在，则块712-722被重复。否则该过程结束。一旦每个处理器已经完成，共享存储器510就将包含整个勘探区域的三维图像。
正如以前提到的，用于给定炮集的处理优选地被分配给单个处理器以此使处理器之间的交叉通信降至最少。然而，对于其中存储需求超过单个处理器的存储能力的任务来说，区域分解可以用来在多个处理器之间分配计算任务。这种分工甚至连同其伴随的对于处理器间通信的需求将优先于由超过单个处理器的存储能力产生的I/O延迟。区域分解技术尤其适用于这里所描述的离散化波动方程，尽管预计会有某些增加的开销(大约为总计算成本的10-15％)。
在某些实施例中，区域分解根据区域来实施。也就是，每个处理器512被分配用于传播和成像的波场体积的一部分。在每个时间步长，每个处理器发送并接收压力场值，用于沿着具有被其它处理器处理的区域的边界的层。该层的厚度是n/2，其是空间差分滤波器的宽度的一半。
下面的计算是为炮集的处理确定存储需求的实例。源或接收器波场的一个时间步长的存储需求是Nb＝NxNyNz(浮点数的大小)(10)其中Nx=L/2+2aΔx+2(Nxp+n2)---(10a)]]>Ny=W/2+2aΔy+2(Nyp+n2)---(10b)]]>
Nz=DΔz+2(Nzp+n2)---(10c)]]>其中L是接收器阵列的长度，W是接收器阵列的宽度，D是将要成像的最大深度，α是用来减小边缘效应的孔径大小，n是空间滤波器的阶，以及Nxp，Nyp，Nzp是x，y和z边界上用于吸收的填充层的厚度。作为一个实例，填充层厚度可从10变至23个格点。对于6km孔径、12km勘探深度和25米空间栅格间距来说，一个时间步长的存储需求可能是大约0.5GB。
因为需要多个时间步长并且正在构建分图像，在这个实例中，用于传播和成像的总存储需求可能是大约2GB。由于源波场的许多时间步长的保存，共享的存储空间明显更大。假定2毫秒时间步长以及10和50之间的抽选因子，每个处理器的总共享存储空间可能是大约60-150GB。
这里所公开的发明提供了在其高度处理磁倾事件、回转波及其复合能力方面对基于单向波动方程的方法的改进。在优选实施例中，本发明利用支持双向波传播同时维持计算可行性的离散波动方程。空间有限差分滤波器的使用还允许区域分解，使得较大的成像工作通过多处理器计算机有效地实施。具有基于照明补偿的成像条件的互相关通过来自多个反射的能量的结合提供了改进的成像。
一旦充分理解了上述公开内容，许多的变化和修改对本领域的技术人员来说将变得显而易见。意图是，下面的权利要求可被解释为包含所有这样的变化和修改。
权利要求
1.一种叠前成像系统，其包含多个处理器，每个处理器被配置成实施源事件的基于双向波动方程的向前波传播以此利用炮集的高阶空间有限差分获取源波场；实施接收波的基于双向波动方程的向后波传播以此利用炮集的高阶空间有限差分获取接收器波场；以及使所述源和接收器波场相关以此获取图像区域。
2.如权利要求1所述的系统，还包含一个或多个共享存储装置，耦合至多个处理器并被共同配置成可存储多个炮集、速度立方体、多个源波场和图像区域。
3.如权利要求2所述的系统，其中按时间方向抽选所述存储的源波场的十分之一。
4.如权利要求1所述的系统，其中利用高阶空间有限差分滤波器来实施所述向前和向后传播以此估计拉普拉斯算子，所述空间滤波器具有比2更高的阶。
5.如权利要求4所述的系统，其中所述空间滤波器具有比6更高的阶。
6.如权利要求4所述的系统，其中所述空间滤波器具有利用级数展开逼近拉普拉斯算子而确定的系数。
7.如权利要求1所述的系统，其中利用时间步长来实施所述向前和向后传播，所述时间步长利用空间栅格间距、最大声速和所述空间滤波器的系数来确定。
8.如权利要求1所述的系统，其中利用时间步长来实施所述向前和向后传播，所述时间步长利用空间栅格间距和最大声速来确定。
9.如权利要求1所述的系统，其中所述多个处理器的每一个还被配置成可使源和接收器波场的边界值与处理相同波场的至少一个其它处理器交换。
10.如权利要求1所述的系统，其中所述多个处理器的每一个还被配置成可通过将所述图像区域与其它处理器确定的任何重叠的图像区域叠加而将所述图像区域贡献给整个勘探图像。
11.一种叠前成像方法，其包含获取用于勘探区域的多个炮集；获取用于所述勘探区域的速度立方体；确定传播时间步长；将每个炮集分配给至少一个处理器，所述至少一个处理器利用离散的高阶空间有限差分波动函数实施源小波的向前传播以此获取用于所述炮集的源波场；利用离散的高阶空间有限差分波动函数实施接收波的向后传播以此获取用于所述炮集的接收器波场；以及使所述源和接收器波场相关以此形成分图像。
12.如权利要求11所述的方法，还包含所述至少一个处理器，使所述分图像与重叠的分图像叠加以此形成整个勘探图像的一部分。
13.如权利要求11所述的方法，其中每个炮集被分配给单个处理器。
14.如权利要求11所述的方法，其中每个炮集被分配给多个处理器，所述多个处理器的每一个在每个时间步长交换多层波场边界值以此使得能够为下一时间步长计算近边界的值。
15.如权利要求11所述的方法，其中利用高阶空间有限差分滤波器来实施所述向前和向后传播以此估计拉普拉斯算子，所述空间滤波器具有比2更高的阶。
16.如权利要求15所述的方法，其中所述空间滤波器具有利用拉普拉斯算子的级数展开逼近而确定的系数。
17.如权利要求16所述的方法，其中利用空间栅格间距、最大声速和所述空间滤波器的系数来确定所述传播时间步长。
18.如权利要求11所述的方法，其中所述组合是所述源和接收器波场的零滞后互相关，其通过所述源波场的能量被标准化。
19.如权利要求11所述的方法，其中所述实施向前传播包括保存所述源波场的均匀时间间隔的型式。
全文摘要
公开了一种用于叠前成像的系统和方法，其包括多个处理器，其中每个处理器被配置成可实施源事件的基于双向波动方程的向前波传播以此利用炮集的高阶空间有限差分获取源波场，并且可实施接收波的基于双向波动方程的向后波传播以此利用炮集的高阶空间有限差分获取接收器波场，以及使源和接收器波场相关以此获取图像区域。
文档编号G01V1/22GK101017204SQ200710084030
公开日2007年8月15日 申请日期2007年2月9日 优先权日2006年2月9日
发明者M·H·卡拉津西尔 申请人:Pgs地球物理公司