CN102395902A - 使用快速面向目标照明计算的地震成像系统及方法 - Google Patents

Abstract

在各种公开的地震成像系统和方法中，使用快速面向目标照明计算技术来获得近似照明值的数据体或“矩阵”。通过这些照明值可以发现“真实反射率”值的图像矩阵。从格林函数获得照明值，格林函数不是以逐个炮点为基础计算或重新计算的，而是以多炮点组进行计算，并与滚动和相结合，从而大大降低计算量。因此，与那些依赖于常规三维波动方程照明的系统相比，所公开的系统和方法能够更加快速地和/或以更好的质量提供目标区域照明。

Description

使用快速面向目标照明计算的地震成像系统及方法

技术领域

地震学用于需要地质信息的勘探、考古研究和工程项目。勘探地震学所提供的数据在与其他可利用的地球物理学、钻井和地质数据联合使用时，能够提供关于岩石类型的结构和分布及其含量的信息。这种信息对于搜寻水、地热储层以及诸如碳氢化合物和矿石的矿床具有极大的帮助。大多数石油公司依靠勘探地震学来选择钻探石油探井的地点。

背景技术

传统的地震学使用人工产生的地震波来勘查地下结构。地震波从源处向下传入地球并从地下结构之间的分界处反射。地面检波器(receiver)探测并记录反射的地震波以供后续分析。尽管直接检验记录的信号可以察觉到一些大尺度结构，但必须对记录的信号进行处理，以消除畸变并且展现地下图像中细微的细节。现有的各种处理方法不能充分消除畸变，并且需要过长的计算时间。本文公开了改进的系统和方法。

附图说明

结合下列附图考虑接下来的详细描述，可以更好地理解各个公开的实施例，其中：

图1显示了例证性的海洋地震勘测环境的侧视图；

图2显示了例证性的海洋地震勘测环境的顶视图；

图3显示了指定通道接收到的交替激发(flip-flop)炮点所产生的例证性的中心点模式(midpoint pattern)；

图4显示了例证性的地震勘测记录系统；

图5显示了例证性的踪迹集；

图6显示了例证性的三维数据体；

图7显示了例证性的炮点几何位置；

图8显示了例证性的地震成像方法的流程图；以及

图9显示了例证性的成像系统。

尽管本发明适用于多种修改及替代形式，但在附图中以示例方式显示了其特定实施例，并在本文中对这些特定实施例进行了详细描述。然而，应理解，附图及详细描述不是要将所公开的实施例限制于所显示的特殊形式，相反，意图是覆盖落入到所附权利要求的实质和范围内的所有修改、等同或替代形式。

具体实施方式

本文公开了使用快速面向目标照明计算技术的多种地震成像系统和方法。计算近似照明值的数据体或“矩阵”，并用其来估计真实反射率值的图像矩阵。从格林函数(Green’s function)获得照明值，格林函数不是以逐个炮点为基础计算或重新计算，而是以多炮点组进行计算，并与滚动和(rolling-sum)相结合，从而大大降低计算量(computationaloverhead)。因此，与那些依赖于常规三维波动方程照明的系统相比，本发明所公开的系统和方法能够更加快速地和/或以更好的质量提供目标区域照明值。

在至少一些实施例中，本发明所公开的地震成像系统包括至少一个存储装置、存储器和至少一个处理器。存储装置存储来自于指定勘测区域的地震勘测的炮点道集(gather)。存储器存储地震成像软件，该地震成像软件在被处理器运行时配置系统，以从存储装置获得多炮点道集并将那些具有重叠的检波器位置的炮点道集分组以形成多炮点道集。对于多炮点道集中的每个炮点和检波器位置，系统同时确定勘测区域中的潜在散射点处的格林函数波场。(如本文进一步详述的，同时获得这些格林函数能够避免大量的重新计算)。系统将格林函数组合，以发现每个潜在散射点处的照明值，并且至少部分基于照明值来计算相应的地震反射率值。系统根据这些反射率值来产生勘测区域的照明，其能够接受进一步的处理或者进行显示以供分析者使用。

在一些实施方式中，每个多炮点道集对应于海洋地震勘测中的一条航海线。通过频域内的小波束分解(beamlet decomposition)和直接局部平面波分解，使用单程波动方程传播来确定格林函数，从而可以得到其他的效果。每次对勘测区域的一个深度切片执行每个格林函数的同时确定，并在格林函数传播到下一深度切片之前执行照明值确定。如上文所述，如果在滚动和的基础上计算所需要的求和项(summationterm)的其中之一，则可以显著节省照明值确定的计算量。通过将照明值的计算限制在感兴趣的目标区域，而不是对整个勘测体积执行计算，可以实现进一步的节省。下面描述所公开的成像系统和方法的上述和其他变形。

通常在陆地上和海上都执行勘探地震。在海上，地震勘测船如图1所示在船的后面配置拖缆。当船向前(沿箭头102的方向)移动时每个拖缆110跟在船100的后面，每个拖缆包括多个均匀间隔的检波器114。每个拖缆110可以进一步包括可编程换向器(diverter)118和可编程深度控制器，其将拖缆从船只航线向外拖到工作偏移距(见图2)并向下拖到所需要的操作深度(图1)。

拖缆110可以长达数千米，通常被构造很多段，每一段的长度为25至100米，包括多达35个检波器或者更加均匀间隔的检波器的组。每个拖缆110包括用于使检波器114和船100上的地震设备互相连接的电缆布线或光缆布线。数据在检波器114附近数字化，并通过电缆以每秒钟7(或更高)百万比特数据的速率传输到船100。

如图1所示，地震勘测船100还可以拖带一个或多个源112。源112可以是脉冲源或振动源。海洋地震学中所使用的检波器114通常被称为水下地震检波器(hydrophone)，通常用压电式换能器构造。可以使用各种合适类型的水下地震检波器，例如盘式水下地震检波器和圆柱形水下地震检波器。源112和检波器114通常配置在海面104以下。船上的处理设备控制源和检波器的操作，并且记录所采集的数据。

地震勘测提供用于在海面104以下成像以展现地下结构(诸如位于海底以下108的结构106)的数据。分析者可以绘制出地下层的形貌并且研究所记录的地震数据的特征，对数据的研究用来确定石油和/或天然气储层的位置。

为了对地下结构106成像，源112发射地震波116，地震波在由于地下结构106(和其他地下结构)而存在声阻抗变化的地方发生反射。通过检波器114的模式来探测反射波。通过记录(作为时间的函数)已从源112传播到地下结构106再到检波器114的到达地震波116，在进行适当的数据处理之后可以获得地下结构106的图像。

图2显示了拖曳一组拖缆110和两个源112的地震勘测船100的俯视图(未按比例)。当船100向前移动时，以所谓的交替激发模式(flip-flop pattern)交替地触发源。使用可编程换向器来提供拖缆之间大致均匀的间隔。拖缆上的指定位置处的检波器与共同的场文件踪迹号(common field file trace number)或共同的通道(common channel)202相关联。

图3显示了两个炮点的例证性的源和检波器的俯视图。对于第一炮点而言，在位置302触发一个源，检波器阵列的所示部分处于位置304(用虚线轮廓线显示)。对于第二炮点而言，在位置306触发源，检波器阵列的所示部分处于位置308(用实线轮廓线显示)。假设目前反映的地下结构是水平的，则到达十二个检波器中的每一个检波器的地震波从源和检波器位置之间的中心点下方的位置反射。因此，第一炮点从十二个中心点308的下方产生反射(用具有垂直交叉阴影线的虚线轮廓线显示)，而第二炮点从十二个中心点310的下方产生反射(用具有水平交叉阴影线的实线轮廓线显示)。作为一个示例，向量312显示地震能量从炮点302到中心点314的传播，相等长度的向量316显示传播到检波器位置的反射地震能量。对于第二炮点306而言，向量318和向量320显示相似的传播路径。应注意，中心点314是被多个炮点“击中”的中心点的其中之一，从而当来自于炮点的信息被处理和组合时，可以获得来自于这些区域的更多的信号能量。地震勘测(用于陆地和海洋两者)通常被设计为提供均匀分布的中心点网格，从而每个中心点均有相当高的平均命中次数。

图4显示了例证性的地震勘测记录系统，其具有连接到总线402以将数字信号传送到勘测船100上的数据记录电路406的检波器114。位置信息及其他参数传感器404也连接到数据记录电路406，以使数据记录电路能够存储对解释所记录的数据有益的附加信息。作为例证，这种附加信息可以包括阵列方向信息和速度信息。

通用数字数据处理系统408被显示为连接到数据记录电路406，进一步被显示为通过总线402连接至定位装置410和地震源112。处理系统408配置记录电路406、定位装置410和地震源112的操作。记录电路406将高速数据流从检波器114采集到诸如光盘或磁盘存储阵列的非易失性存储介质上。定位装置410(包括可编程换向器和深度控制器)控制检波器114和源112的位置。

图4的地震记录系统可包括未在此处具体显示的另外的元件。例如，每个拖揽110可以具有用于连接到数据记录电路的独立总线402。处理系统408包括具有图形显示和键盘或其他接受用户输入的方法的用户接口，处理系统408可以进一步包括网络接口，该网络接口用于将存储的地震勘测数据传送到具有用于处理地震勘测数据的强大计算资源的中央计算设施。

图5显示了检波器114探测和采样的例证性的地震信号(也被称为“踪迹”)。信号指示随时间变化的地震波能量的某种度量(例如位移、速度、加速度、压强)，以可编程采样率在高分辨率(例如24位)下对其进行数字化。可以以不同的方式对这种信号分组，当这样分组时，其被称为“道集”。例如，“共中心点道集”是在定义区域内具有中心点的一组踪迹。“炮点道集”是对地震源的单次引爆(firing)记录的一组踪迹。“多炮点道集”是一组炮点道集，在海洋地震勘测中通常包括沿航海线记录的全部踪迹。

当保持成图5的格式时，所记录的地震勘测数据用处不大。尽管可以用展现大尺度地下结构的曲线图并排标绘制各种记录波形，但是这种结构是变形的，甚至无法看见细微的结构。因此，对数据进行处理，以产生数据体，即数据值的三维阵列，例如如图6所示。数据体表示整个勘测区域内的某种地震属性。三维阵列包括大小均匀的单元，每个单元具有表示该单元的地震属性的数据值。可以表示各种地震属性，在一些实施例中，每个单元具有多个数据值，以表示多个地震属性。合适的地震属性的示例包括反射率、声阻抗、声速和密度。体数据格式本身更容易进行计算分析和视觉渲染，因此，可将数据体称为勘测区域的“三维图像”。

图7显示了各个参数是如何与例证性的炮点的几何位置相关联的。地震能量从地震源沿射线702传播到目标分界面704，并沿射线706向检波器反射。在反射点(在别处用(x，y，z)坐标表示，缩写为向量

)，表面具有法向量

在具有水平面内的x轴和y轴以及朝上的z轴的预定坐标系统中，法向量与每个轴形成角度。法向量

与x轴之间的角度为θ_X，法向量

与y轴之间的角度为θ_Y，法向量

与z轴之间的角度为θ_n。入射射线702和反射射线706相对于法向量成相等(但相反)的“张”角θ_r。入射射线702相对于z轴形成入射角θ_s，反射射线相对于z轴形成接收(又叫做散射)角θ_g。根据下面的关系式，可将入射角和散射角对(θ_s·θ_g)转换为法向角和张角对(θ_n·θ_r)：

${\stackrel{\→}{\mathrm{\θ}}}_n=({\stackrel{\→}{\mathrm{\θ}}}_s+{\stackrel{\→}{\mathrm{\θ}}}_g)/2,$ ${\stackrel{\→}{\mathrm{\θ}}}_r=({\stackrel{\→}{\mathrm{\θ}}}_s-{\stackrel{\→}{\mathrm{\θ}}}_g)/2.---\left(1\right)$

根据炮点的几何位置的这种理解，我们现在开始讨论照明分析及其在地下区域的成像中的作用。宽泛地说，照明分析包括确定有多少地震能量能够被指定体积在检波器的方向上反射。同样地，照明是反射点上的每次击中的炮点几何位置的函数。

照明分析是用于研究图像质量上的采集孔径(acquisition aperture)和重叠结构(overlaying structure)的强大工具。大多数用于预测照明度分布的现有技术是基于射线跟踪模拟。例如参见Bear，G，LU，C.，Lu，R.和Willen，D.，2000，“The construction of subsurface illumination andamplitude maps via ray tracing”，The Leading Edge，19(7)，726-728；Muerdter，D.和Ratcliff，D.，2001，“Understanding subsalt illuminationthrough ray-trace modeling，Part 1：Simple 2D salt models”，The LeadingEdge，20(6)，578-594；Muerdter，D.，Kelly，M.和Ratcliff，D.，2001，“Understanding subsalt illumination through ray-trace modeling，Part 2：Dipping salt bodies，salt peaks，and nonreciprocity of subsalt amplituderesponse”，The Leading Edge，20(7)，688-687；Muerdter，D.和Ratcliff，D.，2001，“Understanding subsalt illumination through ray-trace modeling，Part3：Salt ridges and furrows，and impact of acquisition orientation”，TheLeading Edge，20(8)，803-816；Schneider，W.A.和Winbow，G.A.，1999，“Efficient and accurate modeling of 3-D seismic illumination”，SEGExpanded abstracts 18，633-636。

尽管射线跟踪成本并不高，但是由于所涉及的高频近似以及射线跟踪所产生的奇异性问题，所得到的照明图(illumination map)可能在复杂区域中具有大的误差。参见例如Hoffmann，2001，“Illumination，resolution and image quality of PP-and PS-waves for survey planning”，The Leading Edge，20(9)，1008-1014。为了对付这些问题，使用小波束分解和直接局部平面波分解来得到频域波场的局部角度域信息。更具体而言，使用Wu，R.S.和Chen，L.2002，“Mapping directional illuminationand acquisition-aperture efficacy by beamlet propagator”，SEG ExpandedAbstracts 21，1352；以及Xie，X.B.和Wu，R.S.，2002，“Extracting anglerelated image from migrated wavefield”，Expanded abstracts，SEG 72ndAnnual Meeting，1360-1363所公开的技术来确定作为位置、频率以及入射或散射角的函数的源和检波器波场。这些分解使得能够使用基于波动方程的方法(特别包括单程波动方程传播子)来执行照明分析。关于这种方法的更多细节可在Xie，X.B，Jin，S和Wu，R.S.，2003，“Three-dimensional illumination analysis using wave equation basedpropagator”，SEG Expanded Abstracts 22，989中获得。这种方法的现有实施方式已被证明计算量太大并且是缓慢的，原因在于其逐个炮点来计算照明的方法。

本发明公开的用于计算照明的方法可以使用下列三个特征中的一个或多个特征来加速照明计算：计算目标区域而不是整个数据体的照明；多炮点多炮点地来执行计算，而不是逐个炮点执行计算；在照明计算中使用滚动和，而不是每次都从头开始重新求和。

最终目标是生成地震反射率的数据体，之后作为用来定位储层和矿床的工作的一部分，分析者可以对其进行研究。如Luo，M和Xie，X.B.，2005，“Amplitude recovery from back propagated waves to true scatteredwaves”，Technical Report No.12，Modeling and Imaging Project，University of California，Santa Cruz，25-34所述，“真实反射率”图像可写成：

$R\left(\stackrel{\→}{x}\right)=\underset{{\stackrel{\→}{\mathrm{\θ}}}_s}{\mathrm{\Σ}}\underset{{\stackrel{\→}{\mathrm{\θ}}}_g}{\mathrm{\Σ}}\underset{\mathrm{\ω}}{\mathrm{\Σ}}\underset{s=1}{\overset{n_{\mathrm{stk}}}{\mathrm{\Σ}}}\frac{u_g(\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_g)\·u_s{(\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_s)}^{*}}{n_{\mathrm{stk}}{\left|u_s(\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_s)\right|}^2\underset{r}{\mathrm{\Σ}}\left|G(\stackrel{\→}{x},r,\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_g){\∂}_{n}G(r,\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_g)\right|},---\left(2\right)$

其中

表示反射率(散射系数)，

表示散射点

和频率ω处的从炮点索引号(index number)s传播的具有入射角

的(频域)压强波场，

表示散射点

和频率ω处从检波器索引号r向后传播的具有散射角

的(频域)压强波场。

为炮点s的从散射点

到检波器r的格林函数(频率ω处，散射角为

)，

代表从检波器r向后传播到散射点

的格林函数的角度分量。折叠次数

表示每一点

和入射角为

和接收角为

的频率ω处的倍数(又叫做“多重”、“叠加次数”或“命中次数”)。

在进行特定简化假设以降低式(2)的计算复杂度时，我们首先注意到，可将格林函数看成在时间和空间中响应于脉冲而产生的波场。因此，可以用格林函数和恒定的频率系数来表示源波场，即：

$u_s(\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_s)=f\left(\mathrm{\ω}\right)\·G(\stackrel{\→}{x},\mathrm{\ω},s,{\stackrel{\→}{\mathrm{\θ}}}_s).---\left(3\right)$

这种观察使得可以用格林函数重写式(2)。

通过简化，如果我们忽略1)频率系数f(ω)，2)

和

之间的差，以及3)

和之间的差，则可以用下面的照明矩阵代替式(2)中的分母：

$D(\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)=\underset{s=1}{\overset{n_{\mathrm{stk}}}{\mathrm{\Σ}}}{\left|G(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s)\right|}^2\underset{r}{\mathrm{\Σ}}{\left|G(\stackrel{\→}{x},r,\mathrm{\ω},s,{\mathrm{\θ}}_g)\right|}^2.---\left(4\right)$

为了方便起见，我们将式(2)的分子定义为局部角度域中的卷积图像矩阵：

$I(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)=u_g(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)\·u_s{(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s)}^{*}.---\left(5\right)$

“真实反射率”图像条件变成：

$R(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)=\frac{I(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)}{D(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)}.---\left(6\right)$

因此照明矩阵D(...)充当近似校正因子，以将传统的图像矩阵I(...)转换为真实反射率图像矩阵R(...)。

另外，我们注意到，可以使用式(1)的关系来在整个法向(倾)角θ_n或者整个张角θ_r上求和，以得到倾角索引的照明向量或张角索引照明向量。即，

$D(\stackrel{\→}{x},\mathrm{\ω},{\stackrel{\→}{\mathrm{\θ}}}_n)=\underset{{\mathrm{\θ}}_r}{\mathrm{\Σ}}D(\stackrel{\→}{x},\mathrm{\ω},{\stackrel{\→}{\mathrm{\θ}}}_n,{\mathrm{\θ}}_r),---\left(7\right)$

$D(\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_r)=\underset{{\stackrel{\→}{\mathrm{\θ}}}_n}{\mathrm{\Σ}}D(\stackrel{\→}{x},\mathrm{\ω},{\stackrel{\→}{\mathrm{\θ}}}_n,{\mathrm{\θ}}_r).---\left(8\right)$

现在我们来看图8，其为使用快速面向目标照明计算的例证性的地震成像方法的流程图。在步骤802中，计算机获得地震勘测数据。可如上文所述来采集的地震勘测数据通常以结构化的形式存储在磁盘阵列中，以通过并行处理计算机系统实现高带宽可达性。在步骤804中，计算机获得速度体，即兴趣区域中的每一点的地震速度的指示。可以使用已知的技术从地震勘测数据得到速度体。例如，参见Jon F.Claerbout，Fundamentals of Geophysical Data Processing，p.246-56，其以引用的方式并入本文。或者，速度体可以至少部分得自其他源，例如钻井日志。地震踪迹的初始叠加和分析揭示兴趣区域中的地下结构的大体轮廓，这种大体轮廓可以用于确定需要准确图像的特殊目标区域。

在步骤806中，根据可用的计算资源，将勘测数据分成多炮点组。下面将会进一步说明，使用较小的组需要较少的信息存储空间，但是增加计算负担。如果有足够的存储空间，则选择组的大小，以最大限度地利用重叠的炮点和/或检波器位置。在海洋勘测中，优选的多炮点组是一条或多个航海线。在陆地勘测中，优选的多炮点组覆盖检波器阵列的一个位置的全部炮点，不过这在使用滚动策略的勘测中可以扩展，所述滚动策略是在每次改变一小部分检波器的位置时都采用新的炮点。通过计算机识别的每个多炮点组的特征包括炮点数量、每个炮点的检波器(即踪迹)数量以及源和每个检波器的位置。

在步骤808中，计算机计算从目标区域的每一点到多炮点中的每个炮点和检波器位置的格林函数，即，

$G(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s)$ 和 $G(\stackrel{\→}{x},r,\mathrm{\ω},s,{\mathrm{\θ}}_g),---\left(9\right)$

对目标区域中的所有和多炮点中的所有s和r。如上文所述，通过在源或检波器的位置施加时间脉冲，并且将波动函数传播到目标区域中的每一点

可以确定格林函数。可以使用上文所述的小波束和角度分解方法来进行这种有效的确定。

在步骤810中，对于目标区域中的每一点

计算机使用存储的格林函数来计算式(4)的照明矩阵。但是，我们不是重新每一次都从头开始计算照明矩阵，而是定义了求和项：

$A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)=\underset{r}{\overset{n_{\mathrm{recv}}}{\mathrm{\Σ}}}{\left|G(r,\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)\right|}^2.---\left(10\right)$

在对一个炮点执行该求和之后，通过以滚动的方式“更新”总和，我们就可以利用两个相邻的炮点通常存在许多重叠的检波器的这个事实。这种更新去除了来自于仅在前面的炮点中出现的那些检波器位置的贡献，并且增加了来自于在前面的炮点中不存在但是在当前炮点中存在的那些检波器位置的贡献，即，

$A(\stackrel{\→}{x},\mathrm{\ω},s+1,{\mathrm{\θ}}_g)=A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)-\underset{r}{\overset{n_{\mathrm{sub}}}{\mathrm{\Σ}}}{\left|G(r,\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_g)\right|}^2+\underset{r}{\overset{n_{\mathrm{add}}}{\mathrm{\Σ}}}{\left|G(r,\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_g)\right|}^2.---\left(11\right)$

由于去除和增加的格林函数的数量通常大大小于检波器的总数，因此该“滚动和”方法显著节省计算量。

计算机可以存储所计算的作为位置、频率、炮点索引、入射角和散射角的的函数的照明矩阵值。掌握了照明矩阵，计算机就能够在以后对目标区域中的每一点确定式(6)的反射率条件，从而得到一个炮点对反射率图像的贡献。这些贡献可以在全部炮点上累积，以得到整个目标区域上的累积反射率图像。可以单独地执行这种计算。

在步骤816中，计算机估计是否存在更多的散射角，如果存在，则在步骤818中使散射角θ_g增加，并且循环回到步骤810。以这种方式，步骤810-818实现式(2)中的一个求和。在步骤820中，计算机确定是否存在更多的入射角，如果存在，则在步骤822中使入射角θ_s增加，并且循环回到步骤810。以这种方式，步骤820-822实现式(2)中的另一个求和。类似地，步骤824和步骤826提供在全部炮点索引s上的求和，而步骤828和步骤830提供在整个频率范围ω上的求和。最后，在步骤832中，计算机确定是否存在更多的多炮点，如果存在，则在步骤834中使递增到下一多炮点组。计算机循环回到步骤808，以计算新的多炮点的格林函数。

一旦在步骤832中完成全部循环，则可以计算目标区域中的每一点的总的照明强度并且在步骤836中进行显示。所存储的照明矩阵值还可以为了将来的处理而存储，例如为了确定使用式(6)的近似的反射率图像。

上述方法在三个方面具有突出的优点。第一，前述方法通过使计算的格林函数的数量最小化，充分利用了计算时间和资源。为了进行对比，讨论下列示例，每个示例假设了具有200个航海线的三维勘测，每个航海线具有200个炮点。此外，对每个炮点采集200个检波器踪迹，平均折叠次数(命中次数)大约为100。地下体积使得用于一个格林函数的一个深度切片填充2MB的存储器，而用于格林函数波场的整个数据体填充磁盘上的1GB。假设的时间单位是用于计算一个格林函数数据体的时间。在下表中显示下文所讨论的每个例证性方法的对比结果：

方法	线	炮点/线	检波器	折叠	格林	存储器	磁盘	时间
									1	200	200	200	100	80,000	160G	0	80,000
2	200	200	200	100	80,000	0.002G	80,000G	80,000
									3	1×200	1×200	200	100	200×40,000	0.4G	0	8,000,000
4	2×100	200	200	40	2000×100	4.0G	0	200,000

在第一个例证性方法中，同时计算勘测当中的每个炮点和检波器位置的全部格林函数。也就是说，每个深度层是从前一层传播来的，同时地进行反射率计算，从而不必存储完整的格林函数数据体。该方法提供了高速的优点，但是，从上表中可以观察到，该方法的存储器需求过大。

在第二个例证性方法中，计算并在磁盘中存储勘测当中的每个炮点和检波器位置的全部格林函数。通过这种方法，为每个格林函数保存了完整的波场数据体。该方法和第一个方法一样快速，但是要注意，现在所需要的磁盘容量过高。

在第三个例证性方法中，逐个炮点计算照明矩阵。也就是说，对指定炮点中的炮点和检波器的当前位置计算格林函数，然后对下一炮点重新进行计算。存储器和存储需求非常合适，但是该方法所需要的时间过长。

在第四个例证性方法中，如图8所示来计算照明矩阵，即多炮点多炮点地计算。也就是说，同时计算许多炮点的格林函数波场数据体(或者在该示例中，深度切片波场)，以得到这些炮点的照明。格林函数的数量取决于炮点数量和每个炮点的检波器的数量。同时计算的炮点数量取决于计算机的存储器容量。该方法在存储器/存储需求与计算时间之间做到了非常好的折衷。

图8的方法通过使用式(11)的动求和方法，获得了增加的计算时间。如果检波器位置具有重叠p(通常百分比超过95％)，则计算量将节省大约(2p-1)。

图8的方法通过仅对感兴趣的目标区域计算并存储照明和反射率值来获得第三种计算时间节约。由于目标区域通常只呈现总勘测体积的一小部分，因此目标区域方法去除了许多不必要的计算。当然，这种节约高度取决于目标区域的相对大小和勘测体积。

预期可以用软件的形式实现图8所示的操作，软件可以存储在计算机存储器、长期存储介质和/或便携式信息存储介质中。应注意，图8的例证性方法可以以不同的次序执行，不一定是顺序执行的。例如，并行可对地震数据处理充分有益。在一些处理方法实施例中，可以独立地处理来自于不同勘测区域的数据。在其他实施例中，操作可以是“流水线的”，从而通过按照显示的顺序进行操作来确定个别反射率贡献，即使操作全都是同时发生也是如此。可以在例证性方法中添加额外的操作和/或可以省略所显示的若干操作。

图9显示了用于执行包括使用快速面向目标照明计算的地震成像的地震数据处理的例证性的计算机系统900。个人工作站902通过局域网(LAN)904连接到一台或多台多处理器计算机906，多处理器计算机906则通过LAN连接到一个或多个共用存储单元908。个人工作站902充当到处理系统的用户接口，使得用户可以将勘测数据载入系统、从系统获得和显示图像数据以及配置和监视处理系统的操作。个人工作站902可以采用具有图形显示和键盘的台式计算机的形式，其中图形显示以图形方式显示勘测数据和勘测区域的三维图像，键盘使用户能够移动文件并执行处理软件。

LAN 904提供多处理器计算机906之间以及与个人工作站902的高速通信。LAN 904可以采用以太网的形式。

多处理器计算机906提供并行处理能力，以实现地震踪迹信号到勘测区域图像的适当迅速的转换。每个计算机906包括多个处理器912、分布式存储器914、内部总线916和LAN接口920。每个处理器912对分配的输入数据部分进行操作，以产生地震勘测区域的部分图像。与每个处理器912相关联的是为处理器的使用而存储转换软件和工作数据集的分布式存储器模块914。内部总线916通过接口920提供处理器间的通信以及到LAN网络的通信。通过LAN 904可以提供不同的计算机906中的处理器之间的通信。

共用存储单元908可以是使用用于非易失性数据存储的磁盘介质的大的、独立的信息存储单元。为了提供数据存取速度和可靠性，共用存储单元908可被配置为冗余磁盘阵列。共用存储单元908最初存储来自于地震勘测的速度数据体和炮点道集。照明矩阵值和/或反射率图像体可以存储在共用存储单元908中，以供后续处理。响应于来自工作站902的请求，计算机906可以获得图像体数据并将其提供到工作站，以转换为向用户显示的图形图像。

一旦完全理解上述公开内容，对于本领域技术人员而言，多种变形和修改将是显而易见的。下面的权利要求旨在被理解为包含全部这样的变形和修改。

Claims (20)

1.一种地震成像方法，包括：

获得用于勘测区域的多个炮点道集；

将具有重叠的检波器位置的炮点道集分组，以形成至少一个多炮点道集；

对于指定的多炮点道集中的每个炮点位置和每个检波器位置，同时确定所述勘测区域中的潜在散射点处的格林函数；

组合所述格林函数，以发现每个所述潜在散射点处的照明值；

至少部分基于所述照明值来发现每个所述潜在散射点处的反射率值；以及

向用户显示所述勘测区域的图像，其中所述图像是至少部分基于所述反射率值的。

2.根据权利要求1所述的方法，进一步包括：

对其他多炮点道集重复所述确定操作、组合操作和发现操作；以及

累积从多个多炮点道集确定的反射率值，以形成所述勘测区域的所述图像。

3.根据权利要求2所述的方法，其中在海洋地震勘测中，每个多炮点道集表示一个航海线。

4.根据权利要求1所述的方法，其中通过频域内的小波束分解和直接局部平面波分解，用单程波动方程来确定所述格林函数。

5.根据权利要求4所述的方法，其中每次对所述勘测区域的一个深度切片执行每个格林函数的所述同时确定，用于每个深度切片的确定与所述组合相互交错。

6.根据权利要求4所述的方法，其中在执行所述组合之前，在整个目标区域上执行每个格林函数的所述同时确定。

7.根据权利要求1所述的方法，其中指定散射点

处的所述反射率值能够被表示为频率ω、炮点索引s、入射角θ_s和散射角θ_g处的反射率分量的总和，其中每个所述反射率分量能够被表示为卷积图像值

与照明值

之比，所述卷积图像值

能够被表示为向前传播的源波场与向后传播的接收波场之间的频域乘积。

8.根据权利要求7所述的方法，其中照明值

能够被表示为源贡献与检波器贡献之间

的乘积的总和：

$D(\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)=\underset{s=1}{\overset{n_{\mathrm{stk}}}{\mathrm{\Σ}}}{\left|G(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s)\right|}^{2}A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g),$

其中炮点索引号s从1变化到折叠次数

为从炮点s到散射点

的格林函数，检波器贡献能够被表示为：

$A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)=\underset{r}{\overset{n_{\mathrm{recv}}}{\mathrm{\Σ}}}{\left|G(r,\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)\right|}^2,$

其中

为炮点s的从检波器r到散射点

的格林函数(在频率ω处，具有散射角

)。

9.根据权利要求8所述的方法，其中所述组合包括通过在滚动和的基础上确定检波器贡献

来利用不同炮点的重叠的检波器位置。

10.根据权利要求1所述的方法，其中仅对作为所述勘测区域的子集的目标区域中的潜在散射点执行所述组合。

11.一种地震成像系统，包括：

至少一个存储装置，所述存储装置存储来自于指定勘测区域的地震勘测的炮点道集；

存储器，所述存储器存储地震成像软件；以及

至少一个处理器，所述处理器连接到所述存储器，以执行所述地震成像软件，其中所述软件将所述至少一个处理器配置为：

从所述存储装置检索多个炮点道集；

将具有重叠的检波器位置的炮点道集分组，以形成至少一个多炮点道集；

对于每个多炮点道集中的每个炮点位置和每个检波器位置，同时确定所述勘测区域中的潜在散射点处的格林函数；

组合所述格林函数，以发现每个所述潜在散射点处的照明值；

至少部分基于所述照明值来发现每个所述潜在散射点处的反射率值；以及

生成用于向用户显示的所述勘测区域的图像，其中所述图像是至少部分基于所述反射率值的。

12.根据权利要求11所述的系统，其中所述软件进一步配置所述至少一个处理器，以便累积从多个多炮点道集确定的反射率值，以形成所述勘测区域的所述图像。

13.根据权利要求12所述的系统，其中在海洋地震勘测中，每个多炮点道集表示一个航海线。

14.根据权利要求11所述的系统，其中通过频域内的小波束分解和直接局部平面波分解，用单程波动方程来确定所述格林函数。

15.根据权利要求14所述的系统，其中每次对所述勘测区域的一个深度切片执行每个格林函数的所述同时确定，用于每个深度切片的确定与所述组合相互交错。

16.根据权利要求14所述的系统，其中在执行所述组合之前，在整个目标区域上执行每个格林函数的所述同时确定。

17.根据权利要求11所述的系统，其中指定散射点

处的所述反射率值能够被表示为频率ω、炮点索引s、入射角θ_s和散射角θ_g处的反射率分量的总和，其中每个所述反射率分量能够被表示为卷积图像值

与照明值之比，所述卷积图像值

能够被表示为向前传播的源波场与向后传播的接收波场之间的频域乘积。

18.根据权利要求17所述的系统，其中照明值

能够被表示为源贡献与检波器贡献之间

的乘积的总和：

$D(\stackrel{\→}{x},\mathrm{\ω},{\mathrm{\θ}}_s,{\mathrm{\θ}}_g)=\underset{s=1}{\overset{n_{\mathrm{stk}}}{\mathrm{\Σ}}}{\left|G(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_s)\right|}^{2}A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g),$

其中炮点索引号s从1变化到折叠次数

为从炮点s到散射点的格林函数，检波器贡献能够被表示为：

$A(\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)=\underset{r}{\overset{n_{\mathrm{recv}}}{\mathrm{\Σ}}}{\left|G(r,\stackrel{\→}{x},\mathrm{\ω},s,{\mathrm{\θ}}_g)\right|}^2,$

其中

为炮点s的从检波器r到散射点

的格林函数(在频率ω处，具有散射角

)。

19.根据权利要求18所述的系统，其中作为组合的一部分，所述格林函数发现照明值，所述软件配置所述至少一个处理器以通过在滚动和的基础上确定检波器贡献

来利用不同炮点的重叠的检波器位置。

20.根据权利要求11所述的系统，其中所述软件配置所述至少一个处理器，以便仅对作为所述勘测区域的子集的目标区域中的潜在散射点组合所述格林函数以发现照明值。

Images