地球物理数据的迭代反演方法中的伪迹减少

相关申请的交叉参考

本申请要求2010年5月7日提交的标题为ARTIFACT REDUCTION IN ITERATIVE INVERNSION OF GEOPHYSICAL DATA的美国临时专利申请 61/332,463的权益，其全部内容通过引用合并于此。

技术领域

本发明一般涉及地球物理勘探领域，并且更特别地涉及地球物理数据处 理。特别地，本发明涉及减少数据迭代反演中由反演中做出的近似产生的伪 迹（artifacts）。

背景技术

地球物理反演[1，2]尝试找到最优解释观测数据并满足地质和地球物理约 束的地下性质的模型。存在大量众所周知的地球物理反演方法。这些众所周 知的方法属于两种类别即迭代反演和非迭代反演中的一种。以下是两种类别 中每一种的普遍意义的定义：

·非迭代反演——通过假设一些简单的背景模型并基于输入数据更新该 模型而实现的反演。该方法不使用更新模型作为另一反演步骤的输入。对于 地震数据的情况，这些方法一般被称为成像、偏移、衍射层析成像或博恩反 演。

·迭代反演——涉及地下性质模型的反复改进以便找到满意地解释观测 数据的模型的反演。如果该反演收敛，则最终模型将更好地解释观测数据并 且更接近地近似实际地下性质。迭代反演通常比非迭代反演产生更准确的模 型，但是计算起来昂贵得多。

在地球物理学中普遍采用的两种迭代反演方法是成本函数最优化和级数 法。成本函数最优化包括成本函数S(M)的数值相对于模型M的迭代最小化或 最大化，该成本函数S(M)是计算数据与观测数据之间的错配度的度量（这有 时也被称为目标函数），其中计算数据是通过计算机使用当前地球物理性质 模型和源信号在由给定的地球物理性质模型表示的介质中的物理学支配传播 （physics governing propagation）来仿真的。仿真计算可以通过若干数值方法 中的任何一种来完成，所述数值方法包括但不限于有限差分、有限元或射线 追踪。级数法涉及通过散射方程的迭代级数解进行反演（Weglein[3]）。该解 以级数形式写出，其中级数中的每项对应于散射的较高阶。在此情况下的迭 代对应于向该解添加级数中较高阶的项。

成本函数最优化方法是局部的或全局的[4]。全局方法仅仅涉及计算模型 群体{M₁，M₂，M₃，…}的成本函数S(M)以及从该群体中选择近似最小化S(M) 的一组一个或更多个模型。如果希望进一步的改进，则该新选择的模型组可 以被用作基础来生成能够再次相对于成本函数S(M)被再次测试的新模型群 体。对于全局方法，测试群体中的每个模型都可以视为是一次迭代，或者在 更高水平下，每组被测试的群体可以视为是一次迭代。众所周知的全局反演 方法包括蒙特卡罗法、仿真退火（simulated annealing）法、遗传和进化算法。

局部成本函数最优化包括：

1.选择起始模型；

2.计算成本函数S(M)关于描述该模型的参数的斜率；

3.搜索在梯度方向上作为起始模型的扰动的更新模型，该更新模型更好 地解释观测数据。

该进程通过使用新的更新模型作为另一梯度搜索的起始模型来重复。该 过程继续进行直到找到满意地解释观测数据的更新模型。普遍使用的局部成 本函数反演方法包括梯度搜索法、共轭梯度法和牛顿法。

如上所述，因为迭代反演产生更准确的地下参数模型，所以迭代反演比 非迭代反演更优选。不利的是，迭代反演在计算方面是非常昂贵的，以至于 将其应用于许多感兴趣问题是不实际的。这一高计算成本是所有反演技术都 需要许多计算密集的正向和/或反向仿真的结果。正向仿真意味着在时间上向 前的数据计算，而反向仿真意味着在时间上向后的数据计算。

由于其高昂的计算成本，迭代反演经常需要应用加速计算的一些类型近 似。不利的是，这些近似通常在最终反演模型中导致误差，该误差可以作为 在反演中采用的近似的伪迹（artifacts）观察到。

需要的是一种迭代反演数据的一般方法，该方法允许应用近似而不在最 终的反演模型中生成伪迹。本发明满足该需要。

发明内容

物理性质模型给出一个或更多个地下性质作为区域中的位置的函数。地 震波速是一种此类物理性质，而且（例如）密度、p-波波速、横波波速、若干 各向异性参数、衰减（q）参数、孔隙度、渗透度和电阻率也是此类物理性质。 本发明是用于减少地下物理性质模型中的伪迹的方法，除了信源编码之外， 所述伪迹还由迭代的计算机化的地球物理数据反演过程中的近似导致，所述 方法包含随着所述迭代进行而改变所述近似。在一个特别实施例中，本发明 是用于反演测量的地球物理数据以确定地下区域的物理性质模型的计算机实 施的方法，其包含：

（a）假设地下区域的物理性质模型，所述模型提供在遍及所述地下区域 的多个位置处的至少一个物理性质的数值；

（b）选择具有其中计算由对所述物理性质模型的更新组成的步骤的迭代 数据反演过程，所述更新使得所述物理性质模型更符合从所测量的地球物理 数据；

（c）在所述计算中做出近似，所述近似除通过信源编码之外加速所选择 的迭代数据反演过程或者产生准确度折衷；

（d）使用计算机通过所述近似并使用所述物理性质模型执行所选择的迭 代数据反演过程的一个循环；

（e）使用所述计算机执行下一个迭代反演循环，其中做出改变或不改变 所述近似的一些方面的选择；

（f）如必需则重复（e），在一些或全部迭代循环中改变所述近似，直到 满足所选择的收敛标准或达到另一停止条件的最终迭代；以及

（g）从所述最终迭代下载更新的物理性质模型，或将所述更新的物理性 质模型保存到计算机存储装置。

在本发明的一些实施例中，一个或更多个伪迹类型在反演结果中被识别 为由所述近似导致，并且在一些或全部迭代循环中改变的所述近似的方面被 选择为对一个或更多个所识别的伪迹类型的伪迹有影响。对伪迹的影响可以 使得来自一个近似的伪迹不与来自使用具有已改变方面的近似的另一迭代循 环的伪迹建设性地累加。

附图说明

本发明及其优点将通过参考以下详细描述和附图更好地理解，其中：

图1是示出在此公开的一般方法中的基本步骤的流程图；

图2是示出图1的方法的特别实施例中的基本步骤的流程图，其中通过 对信源进行编码和求和来近似目标函数；

图3-5属于图2的发明实施例的示例性应用：

图3示出地震波速模型，从其计算出针对该示例的地震数据；

图4示出使用在图2中概述的反演方法反演来自图3中的地震波速模型 的数据；

图5示出在没有在多次迭代之间改变用于编码信源的代码的步骤的情况 下，使用在图2中概述的反演方法反演来自图3中地震波速模型的数据；

图6是示出图1的方法的特别实施例中的基本步骤的流程图，其中近似 是改变用于数值反演的栅格单元的大小，以便仅在需要时使用精细栅格；

图7-9属于图6的发明实施例的示例性应用：

图7是地震波速模型，从其计算出针对该示例的地震数据；

图8示出使用在图6中概述的反演方法反演来自图7中的地震波速模型 的数据；

图9是在没有在多次迭代之间改变人工反射生成器的深度的步骤的情况 下，使用在图6中概述的反演方法反演来自图7中地震波速模型的数据；

图10是示出图1的方法的特别实施例中的基本步骤的流程图，其中近似 仅使用测量数据的子集；

图11-13属于图10的发明实施例的示例性应用：

图11示出地震波速模型，从其计算出针对该示例的地震数据；

图12示出使用在图10中概述的反演方法反演来自图11中的地震波速模 型的数据；以及

图13示出在没有在多次迭代之间随机改变测量数据的子集的步骤的情况 下，使用在图10中概述的反演方法反演来自图11中的地震波速模型的数据。

由于专利约束，图3-5、7-9和11-13是彩色显示的灰度转换图。

本发明将结合其优选实施例进行描述。然而，就以下详细描述专门针对 本发明的特别实施例或特别用途来说，希望其仅是说明性的，并且不解释为 限制本发明的保护范围。相反，其意图覆盖可以包括在由随附权利要求定义 的本发明的保护范围内的全部替换、修改和等价物。

具体实施方式

本发明是用于在地球物理数据的迭代反演期间减少由近似的应用导致的 伪迹。地球物理反演尝试找到最优解释观测到的地球物理数据的地下性质的 模型。在阐述本发明的方法的整个过程中始终使用地震数据的示例，但该方 法可以有利地应用于地球物理勘探的任何方法和任何类型的地球物理数据。 使用迭代方法最准确地执行数据反演。不利的是，迭代反演在计算上经常是 禁止性地昂贵的。执行地球物理数据的昂贵的正向和/或反向仿真花费了迭代 反演中的大部分计算时间（在此正向的意思是在时间上向前，并且反向的意 思是在时间上向后）。这些仿真的高成本一部分是由于输入数据中的每个地 球物理源必须在仿真软件的独立计算机运行中计算。因此，仿真的成本与地 球物理数据中的源数目成比例，对于地球物理勘查通常为约1000到10000个 源。在通常的实践中，在反演期间应用近似以减少反演的成本。这些近似导 致反演模型中的误差或伪迹。本发明通过在反演的多次迭代之间改变近似的 一些方面以使得在一次迭代期间的伪迹不与其它迭代中的伪迹建设性地累加 来减轻这些伪迹。因此减少了反演模型中的伪迹。

在迭代反演期间做出的导致伪迹的一些普遍近似包括：

1.应用到测量数据的处理

2.仿真中的不准确边界条件

3.仿真中的近似（例如，仿真器中使用的导数的低阶近似或计算中使用 的栅格单元的大小）

4.模型的参数化中的近似（例如，使用过于粗糙而不能准确表示模型中 的变化的参数的空间栅格）。

在地球物理学中普遍采用的两种迭代方法是成本函数最优化和级数法。 本发明可以应用于这两种方法。接下来是这些方法中的每一个的概述。

迭代成本函数最优化

成本函数最优化通过相对于地下模型M最小化成本函数S(M)（有时也被 称为目标函数）的数值来实现，该成本函数S(M)是观测（测量）的地球物理 数据与通过假设模型的仿真计算的相应数据之间的错配度的度量。经常在地 球物理反演中使用的简单成本函数S是：

$S\left(M\right)=\underset{g=1}{\overset{N_g}{\Sigma }}\underset{r=1}{\overset{N_r}{\Sigma }}\underset{t=1}{\overset{N_t}{\Sigma }}{|{\psi }_{\mathrm{calc}}(M,g,r,t,w_g)-{\psi }_{\mathrm{obs}}(g,r,t,w_g)|}^N---\left(1\right)$

其中

N=成本函数的范数（通常使用最小二乘法或L2-范数，在此情况下N= 2），

M=地下模型，

g=集合（gather）指数（对于点源数据，这可以对应于各个源），

N_g=集合的数目，

r=集合内的接收器指数，

N_r=集合中接收器的数目，

t＝数据记录内的时间样本指数，

N_t=时间样本的数目，

ψ_calc=从模型M计算的地球物理数据，

ψ_obs=测量的地球物理数据，以及

W_g=集合g的源特征，即没有大地滤波效应的源信号

方程1中的集合可以是能够在正向建模程序的一次运行中仿真的任何类 型集合。对于地震数据，该集合对应于地震爆破，尽管爆破可能比点源更一 般[5]。对于点源，集合指数g对应于各个点源的位置。对于平面波源，g可 以对应于不同的平面波传播方向。该广义源数据ψ_obs可以在现场获取或者可 以从使用点源获取的数据合成。另一方面，计算数据ψ_calc通常可以在正向建 模时通过使用广义源函数直接计算（例如，对于地震数据，正向建模通常意 味着各向异性粘-弹性波传播方程的解或其一些近似）。对于许多类型的正向 建模（包括有限差分建模），广义源需要的计算时间约等于点源需要的计算 时间。模型M是地下区域的一个或更多个物理性质的模型。地震波速是一种 此类物理性质，而且（例如）p-波波速、横波波速、若干各向异性参数、衰减 （q）参数、孔隙度和渗透度也是此类物理性质。模型M可以表示单个物理性 质，或者它可以含有许多不同的参数，这取决于反演的复杂程度。通常，地 下区域被细分成离散单元，每个单元由每个参数的单个数值表征。

迭代反演的一个主要问题是计算ψ_calc占用大量计算时间，并且因此成本 函数S的计算非常耗时。此外，在典型的反演方案中，必须为许多不同的模 型M计算该成本函数。

迭代级数反演

除成本函数最优化之外，地球物理反演也可以使用迭代级数法来实施。 用于实现该目的的普遍方法是迭代李普曼-施温格（Lippmann-Schwinger）方 程[3]。李普曼-施温格方程描述介质中的波散射，其由感兴趣的物理性质模型 表示为较简单模型的扰动。该方程是用于确定来自感兴趣模型的波散射的级 数展开式的基底，优点是级数仅需要在较简单模型中执行计算。该级数也可 以被反演以从测量数据形成允许确定感兴趣模型的迭代级数，这同样仅需要 在较简单模型中执行计算。李普曼-施温格方程是可以应用于包括地震波的所 有类型地球物理数据和模型的一般形式。该方法开始于两个方程：

LG=-I      （2）

L₀G₀=-I    （3）

其中L、L₀分别是实际和参考差分算子，G和G₀分别是实际和参考格林算子， 并且I是单位算子。注意G是测量的点源数据，并且G₀是来自初始模型的仿 真点源数据。用于散射理论的李普曼-施温格方程是：

G=G₀+G₀VG    （4）

其中V=L–L₀，可以从其求取真实模型与初始模型之间的差值。

通过以下方式针对V迭代地求解方程4：首先将方程4以级数展开（为G 的第一近似假设G=G₀，以此类推），从而得到：

G=G₀+G₀VG₀+G₀VG₀VG₀+…（5）

然后V展开为级数：

V=V⁽¹⁾+V⁽²⁾+V⁽³⁾+…   （6）

其中V⁽ⁿ⁾是作为数据的残数中的第n阶的V的部分（此处数据的残数是在表 面测量的G-G₀）。将方程6带入方程5并合并相同阶的项，得到针对前3 阶的以下方程组：

G-G₀=G₀V⁽¹⁾G₀            （7）

0=G₀V⁽²⁾G₀+G₀V⁽¹⁾G₀V⁽¹⁾G₀（8）

0=G₀V⁽³⁾G₀+G₀V⁽¹⁾G₀V⁽²⁾G₀+G₀V⁽²⁾G₀V⁽¹⁾G₀+G₀V⁽¹⁾G₀V⁽¹⁾G₀V⁽¹⁾G₀（9）

并且对于V中的更高阶有相似的结果。可以通过以下方式迭代地求解这些方 程：首先通过在V⁽¹⁾的两侧上反演G₀来针对V⁽¹⁾求解方程7，从而产生：

V⁽¹⁾=G₀^-1（G-G₀）G₀^-1    （10）

然后将来自方程10的V⁽¹⁾带入方程8，并且针对V⁽²⁾求解该方程从而产生：

V⁽²⁾=-G₀^-1G₀V⁽¹⁾G₀V⁽¹⁾G₀G₀^-1      （11）

对于V的更高阶同样如此。

方程10包括在源和频率上的和，其可以明确写作：

$V^{\left(1\right)}=\underset{\omega }{\Sigma }\underset{s}{\Sigma }{G}_0^{-1}(G_s-G_{0s})G_{0s}^{-1}$

其中G_s是针对源s的测量数据，G_0s是针对源s通过参考模型仿真的数据，并 且G_0s^-1可以解释为来自源s的向下外推的源特征。方程10在频域中实施时可 以被解释如下：（1）通过参考模型向下外推每个源的源特征（G_0s^-1项），（2） 针对每个源，通过参考模型向下外推残余数据的接收器（G₀^-1（G_s-G_0s）项）， （3）将这两个域相乘，然后在全部源和频率上求和。该方法中的向下外推可 以使用地球物理仿真软件例如使用有限差分来执行。

示例性实施例

图1的流程图示出本发明的方法的一个实施例中的基本步骤。在步骤110 中，选择将改进反演过程的一些方面的近似。该改进可能表现为加速的形式 而不是提高的准确度。此类近似的示例包括使用近似目标函数或在仿真软件 中使用近似。经常选择这些近似来减少反演的计算成本。然而，除计算加速 之外，该改进可以替代地产生准确度折衷，即接受计算的一个方面的更高的 不准确度，作为回报得到一些其它方面的更高的准确度。在步骤140中，基 于测量数据130生成对假设的物理性质模型120的更新。在步骤140中，在 110中选择的近似被用于执行更新计算。使用迭代局部成本函数最优化作为迭 代反演的示例，作为该术语在此使用的“更新计算”包括但不限于计算目标 （成本）函数、目标函数梯度以及实现前述各项所需要的所有正向建模。步 骤140产生更新的物理性质模型150，其应该比假设的物理性质模型120的物 理性质更接近实际地下性质。常规地，通过将该更新的物理模型150和测量 数据130反馈回到步骤140中的更新方法来产生进一步改进的物理性质模型， 从而进一步改进该更新的物理性质模型150。该常规迭代反演方法的缺点在于 由在步骤110中选择的近似导致的反演中的任何伪迹可能在反演中建设性地 增强并损害最终反演结果。

不同于直接返回步骤140，本发明的方法插入步骤160，其中在步骤110 中选择的近似的一些方面被改变以使得由该近似导致的伪迹将改变并因此不 被步骤140的迭代增强。通过该措施，将减轻由在步骤110中选择的近似导 致的伪迹。

近似和相应伪迹的示例

以下表格含有步骤110的示例，即可以有利地用在数据反演中并且适于 （步骤160）本发明的应用的近似的示例。该表格的第一列列出可以与本发明 一起使用的近似。第二列列出与每种近似关联的伪迹。最后一列列出可以在 多次迭代之间变化从而导致多次迭代之间的伪迹改变的特征，该改变导致伪 迹不相干地累加到最终反演模型并因此被减轻。

上面的列表不是穷举的。该列表仅包括减少计算时间的近似的示例。 有时在一个区域中折衷不准确度从而在另一区域中获得更高准确度是有 利的。这种准确度折衷类型的近似的示例是在正向建模中使用较不准确的 吸收边界条件以便使得梯度计算更准确。需要吸收边界条件以求解支配波 传播的差分方程，例如在地震数据情况下的各向异性粘-弹性波传播方程 （或其一些近似），或者在电磁数据情况下的麦克斯韦方程。通常，准确 度折衷涉及在方法的一个方面牺牲准确度，并以另一方面提高的准确度作 为回报。

测试示例1——编码目标函数

图2-5表示使用对目标函数的近似执行反演的综合示例，在该近似中 对测量数据中的地震源进行编码然后求和；参见Jerome Krebs等人的美国 专利公开US 2010-0018718。该近似加速了反演，因为可以使用仿真软件 的一次运行来估计被编码的目标函数，而不是如同在常规反演的情况下为 每个源运行一次。图2是将图1集中在该特别实施例上的自解释流程图， 其中步骤210示出编码近似。

该示例中的地球物理性质模型仅是声波波速的模型。图3示出该示例 的基础速度模型，即被反演并用来生成有待反演的数据的“未知”模型。 阴影表示在每个深度和横向位置的速度，如由向右的“彩色”条所示。图 4示出由应用如图2中的流程图概述的本发明所导致的反演。在该示例中， 通过将这些源随机地乘以+1或-1来将它们编码。在步骤260中，通过改变 用来生成代码的随机数种子来改变这些源的编码，其中所述代码用来编码 这些源。注意在图3中示出对基础模型的良好匹配。

图5示出应用在图2的流程图中概述的反演方法的结果，但消除了步 骤260的发明性特征。注意图5中的反演由串扰噪声支配（反演的斑点状 外观），而该串扰噪声伪迹在源自本发明的反演中基本上不可见（图4）。

同时源的编码与将编码从一个迭代变化到下一个迭代的技术一起在 Jerome Krebs等人的美国专利公开US 2010-0018718中公开；参见该专利 公开中的62段和权利要求3。然而，Jerome Krebs等人的美国专利公开 US 2010-0018718没有认识到或公开该编码发明是在此公开的一般发明的 特定示例。

测试示例2——生成人工反射的近似

图6-9图解说明使用对生成人工反射的仿真器的近似来执行反演的综 合示例。此类近似的示例是使用有限差分仿真器以使栅格中单元的大小随 着距离表面的深度而改变。这种近似加速了反演，因为仿真器中的栅格可 以被调整从而以深度变化方式将其最优化。通常，有限差分仿真器的浅部 分所需要的栅格单元小于模型中较深处所需要的栅格单元。由该近似生成 的伪迹是在栅格单元大小的改变之间的边界处的人工反射。

图6是在该示例中图解说明的本发明的实施例的流程图。在该示例中， 可变栅格仿真器实际上不被用来生成人工反射器。相反（步骤610）通过 在500米深度在密度模型中设置虚拟不连续性来生成人工反射。该不连续 密度模型被仿真器用于模型更新，但恒定密度模型被用于生成测量数据 （图6中的630）。然后以使得仅更新速度模型的方式执行反演（640）， 从而贯穿整个反演的迭代保持虚拟密度不连续性。

该示例中的地球物理性质模型仅是声波波速的模型（620）。图7示 出该示例的基础速度模型（被反演并用来生成有待反演的数据的模型）。 阴影表示在每个深度的速度。图8示出由应用如图6的流程图概述的本发 明所导致的反演。在该示例中，在步骤660中，使用以500米为中心并具 有100米的方差的正态分布来随机改变虚拟密度差的深度。注意图7中示 出对基础模型的良好匹配。在图7-9和图11-13中，速度被绘制为等于反 演速度除以初始速度的无量纲相对速度（v/v(init)），该初始速度是速度模 型预期的起始猜想。

图9示出应用图6的流程图中概述的反演方法的结果，但消除了步骤 660的发明性特征。可以注意到图9中的反演在500米深度处具有清晰可 见的人工反射910，而该人工反射在使用本发明的方法的反演中基本上不 可见（图8）。

测试示例3——测量数据的随机子集

图10-13表示使用对测量数据的近似来执行反演的综合示例。此类近 似的示例是使用测量数据的子集（图10中的1010）。该近似减少了测量 数据的量，由此加速反演，因为反演的计算时间与测量数据的数量成正比。 在通常的反演中，需要全部测量数据以维持高的水平分辨率，并因此在通 常实践中不使用该近似。由该近似生成的伪迹是由稀疏源位置和水平分辨 率恶化导致的反演模型中的覆盖区（footprint）。图10是将图1的步骤集 中在该示例中使用的本发明实施例上的流程图。在该示例中，测量数据的 子集（图10中的1030）在该反演中被使用，例如50个测量数据中的5个 数据的子集。

该示例中的地球物理性质模型仅是声波波速的模型。图11示出该示 例的基础速度模型（被反演并用来生成有待反演的数据的模型）。阴影表 示在每个深度的速度。图12示出由应用如图10的流程图概述的本发明所 导致的反演。在该示例中，在步骤1060中，随着反演迭代增加而随机选 择测量数据的子集。这导致数据的不同子集在每个迭代循环中使用。图12 使用百分之十的测量数据示出对图11中示出的基础模型的良好匹配。

图13示出应用图6中的流程图概述的反演方法的结果，但消除了发 明性的伪迹减少步骤1060。可以注意到图13中的反演在2000米以下的较 深部分具有人工覆盖区，并在整个反演模型中具有短波长噪声，而该覆盖 区噪声在使用本发明方法的反演中减轻（图12），并且短波长噪声是不可 见的。

应理解图2、6和10的流程图表示在图1中更一般地描述的本发明的 特定实施例的示例。

前述专利申请针对本发明的特别实施例以便图示说明本发明。然而， 对本领域技术人员明显的是对在此描述的实施例的许多修改和变化是可 能的。希望全部此类修改和变化在由随附权利要求定义的本发明的保护范 围内。本领域技术人员容易认识到在本发明的实际应用中，本发明方法中 的至少一些步骤（通常是步骤140-160，并且经常在120中生成模型）在 计算机上执行，即本发明是计算机实施的。在这种情况下，最终的地下更 新的物理性质模型可以下载或保存到计算机存储装置中。

参考文献

1.Tarantola，A.，“Inversion of seismic reflection data in the acoustic approximation”，Geophysics 49，1259-1266（1984）。

2.Sirgue，L.和Pratt G.“Efficient waveform inversion and imaging：A strategy for selecting temporal frequencies”，Geophysics69，231-248（2004）。

3.Weglein，A.B.，Araujo，F.V.，Carvalho，P.M.，Stolt，R.H.， Matson，K.H.，Coates，R.T.，Corrigan，D.，Foster，D.J.，Shaw，S.A.， 以及Zhang，H.，“Inverse scattering series and seismic exploration”，Inverse Problems 19，R27-R83（2003）。

4.Fallat，M.R.，Dosso，S.E.，“Geoacoustic inversion via local，global， and hybrid algorithms”，Journal of the Acoustical Society of America 105， 3219-3230（1999）。

5.Berkhout，A.J.，“Areal shot record technology”，Journal of Seismic Exploration l，251-264（1992）。

6.Krebs，Jerome等人，“Iterative Inversion of Data from Simultaneous Geophysical Sources”，美国专利申请公开US 2010-0018718（01-28-2010）。

7.Clapp，R.G.，“Reverse time migration with random boundaries”，SEG International Exposition and Meeting（Houston），Expanded Abstracts，2809-2813 （2009）。