利用同时和顺序源方法进行全波形反演的混合方法

摘要

源（或接收器）编码地球物理数据集合的同时全波场反演方法，以确定地下区域的物理性质模型，特别适于在数据采集中不满足固定接收器几何条件的考察。首先，通过同时编码（203）源反演（205）来反演满足固定接收器条件的数据浅层时窗（202）。然后，利用来自浅层时窗（206）的物理性质模型作为起始模型（207），通过稀疏顺序源反演（209）来反演数据的较深时窗（208）。可选地，利用浅层时窗模型模拟漏测远炮检距数据（211），生成满足固定接收器假设的数据组，然后将该数据组进行源编码（212）和通过同时源反演（214）进行反演。

说明书

相关申请的交叉参考

本申请要求2010年9月27日提交的美国临时申请61/386,828的权益，本文 将其被赋予的全部权限引入作为参考。本申请涉及名称为“Simultaneous Source  Encoding and Source Separation as a Practical Solution for Full Wavefield Inversion”并 要求2010年9月27日提交的美国临时申请61/386,831和2010年10月13日提交 的美国专利申请12/903,749的权益的具体申请。本文将这些所涉及的申请被赋予 的全部权限也引入作为参考。

发明领域

本发明总体涉及地球物理勘探领域，更具体地涉及地球物理数据处理。具体 地，本发明是得自多种地球物理源如震源的数据反演方法，包括地球物理模拟—— 在进行一次模拟中计算来自多个同时活动的地球物理源的数据。

发明背景

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

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

迭代反演——涉及重复改进地下性质模型从而建立理想地解释观测数据的模 型的反演。如果反演收敛，则最终模型将更好地解释观测数据并将更紧密地接近 实际的地下性质。迭代反演通常生成比非迭代反演更准确的模型，然而计算起来 昂贵得多。

迭代反演总体上相对于非迭代反演是优选的，因为其生成更准确的地下参数 模型。不幸地，迭代反演的计算费用如此昂贵，以致于将其应用于多种感兴趣的 问题是不切实际的。这种高计算成本是源于所有反演技术均需要多个计算密集型 模拟的事实。任何单独模拟的计算时间与所要反演的源数成比例，并且地球物理 数据一般存在大量源，在此前述使用的术语源是指源设备的激活位置。迭代反演 的问题加剧，因为必须计算的模拟数与反演的迭代数成比例，并且所需迭代数一 般属于成百上千的级别。

用于地球物理学的最常用迭代反演方法是成本函数优化。成本函数优化包括 成本函数S(M)的数值相对于模型M的迭代最小化或最大化，该成本函数S(M)是计 算数据与观测数据之间错配度的度量（其有时也被称为目标函数），其中计算数据 是通过计算机使用当前地球物理性质模型和源信号在由给定的地球物理性质模型 表示的介质中的物理学支配传播（physics governing propagation）来模拟的。模拟 计算可通过几种数值方法中任一种完成，该数值方法包括但不限于有限差分、有 限元或射线追踪。模拟计算可在频率域或时间域中进行。

成本函数优化方法是局部或全局的[3]。全局方法简单地包括计算模型群体 {M₁，M₂，M₃,…}的成本函数S（M）和从大约最小化S(M)的群体选择具有一个或 多个模型的组。如需进一步改进，则此新选定的模型组可被用作生成新模型群体 的基础，该新模型群体可关于成本函数S(M)再次进行测试。对于全局方法，测试 群体中的各模型可被认为是迭代，或在较高水平上各组所测群体可被认为是迭代。 公知的全局反演方法包括蒙特卡罗（Monte Carlo）、模拟退火、遗传和进化算法。

不幸地，全局优化方法一般收敛极其缓慢，因此大多数地球物理反演是基于 局部成本函数优化。算法1概括了局部成本函数优化。

算法1–进行局部成本函数优化的算法。

此程序通过利用新更新的模型作为另一梯度搜索的起始模型进行迭代。该过 程持续，直到发现理想地解释观测数据的更新模型。常用的局部成本函数反演方 法包括梯度搜索、共轭梯度和牛顿法。

在声学近似中地震数据的局部成本函数优化是普遍的地球物理反演作业，并 通常是其他类型的地球物理反演的示例。当在声学近似中反演地震数据时，成本 函数可写作：

$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 }}W({\psi }_{\mathrm{calc}}(M,r,t,w_g)-{\psi }_{\mathrm{obs}}(r,t,w_g))$    （方程式1）

在此：

S=成本函数，

M=N个描述地下模型的参数（m₁，m₂，…m_N）的向量，

g=集合（gather）指数，

w_g=集合g的源函数，其是空间坐标和时间的函数，对于点源而言其是空间 坐标的δ函数，

N_g=集合数，

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

N_r=集合中的接收器数，

t=道（trace）中的时间样本指数，

N_t=时间样本数，

W=最小化标准函数（我们通常选择W（x）=x²，其是最小二乘方（L2）标 准），

ψ_calc=从模型M计算的地震压力数据，

ψ_obs=测量地震压力数据。

该集合可以是可在一轮地震正向建模程序中模拟的任何类型的集合。通常集 合相应于地震炮点（shot），尽管炮点可比点源更普遍（general）。关于点源，集合 指数g相应于各个点源的位置。关于平面波源，g将相应于不同的平面波传播方向。 这种广义源（generalized source）数据ψ_obs可以在现场得到或可由利用点源得到的 数据合成。另一方面，计算数据ψ_calc通常可通过利用正向建模时广义化的源函数 直接计算。关于多种类型的正向建模——包括有限差分建模，广义源所需的计算 时间粗略地等同于点源所需的计算时间。

方程式（1）可简化为：

$S\left(M\right)=\underset{g=1}{\overset{N_g}{\Sigma }}W\left(\delta (M,w_g)\right)$    （方程式2）

其中，现暗指基于接收器和时间样本的总和，并且

$\delta (M,w_g)={\psi }_{\mathrm{calc}}(M,w_g)-{\psi }_{\mathrm{obs}}\left(w_g\right)$    （方程式3）

反演试图更新模型M，以使S(M)是最小的。其可通过局部成本函数优化实现， 该局部成本函数优化如下更新给定模型M^(k)：

$M^{(k+1)}=M^{\left(k\right)}-{\alpha }^{\left(k\right)}{▿}_{M}S\left(M\right)$    （方程式4）

其中，k是迭代数，α是模型更新的标量大小，(M)是关于模型参数取定 的错配度函数的梯度。模型扰动或模型更新值通过目标函数梯度乘以步长α计算， 其必须被反复计算。

基于方程式（2），可得到下列成本函数梯度的方程式：

${▿}_{M}S\left(M\right)=\underset{g=1}{\overset{N_g}{\Sigma }}{▿}_{M}W\left(\delta (M,w_g)\right).$    （方程式5）

因此，为计算成本函数的梯度，必须分别计算各集合对成本函数的贡献梯度， 然后加和那些贡献。因此，计算所需的计算量（computational effort）是 确定单个集合对梯度的贡献所需的计算量的N_g倍。对于地球物理学问题，N_g通常 相应于地球物理源数，并且属于10,000至100,000的级别，大幅扩大了计算的成本。

注意，的计算需要计算关于N个模型参数中的每一个m_i的导数W （δ）。由于对于地球物理学问题而言N通常很大（通常大于一百万），这种计算可 能极其耗时——如果其必须对于每一个单独的模型参数进行。幸运地，可利用伴 随（adjoint）法同时（at once）对所有模型参数有效地进行这种计算[1]。通过如下 算法总结用于最小二乘方目标函数和网格模型参数化的伴随法：

算法2–利用伴随法计算网格模型的最小二乘方成本函数梯度的算法。

虽然利用伴随法计算梯度相对于其他方法有效，但是其仍十分昂贵。具体地， 伴随法需要两种模拟，一种在时间上正向，一种在时间上反向，而对于地球物理 学问题这些模拟通常是极其计算密集型的。而且，如上所述，这种伴随法计算必 须分别对每一个测量数据集合进行，增加了N_g倍的计算成本。

所有反演种类的计算成本可通过自源组合反演数据而非分别反演源而降低。 这可被称为同时源反演。几种类型的源组合是已知的，包括：相干（coherently） 加和近距离的源以生成有效源，该有效源生成一些所需形状的波前（例如，平面 波）、加和宽距离的源、或完全或部分地叠加反演前的数据。

通过反演组合源实现的计算成本降低由于组合数据的反演通常产生不太准确 的反演模型的事实而至少部分抵消。这种准确性的损失是缘于如下事实：当加和 各个源时漏测信息，因此加和数据限制反演模型的强度不如未加和的数据。加和 过程中的这种信息漏测通过在加和前编码各炮点记录被最小化。组合前进行编码 显著保留同时源数据中的更多信息，因此更好地限制反演[4]。编码还允许组合近 间距的源，因此允许给定计算区域的更多源组合。多种编码方案可用于这种技术， 包括时移编码和随机相位编码。此背景部分的其余内容简要地评述多种已公开的 地球物理同时源技术——编码和非编码的。

Van Manen [6]提出利用地震干涉测量方法加速正向模拟。地震干涉测量通过 将源布置在目标区域的边界各处而工作。对这些源分别建模，并且记录需要格林 函数（Green’s function）的所有位置的波场。任两个记录位置之间的格林函数然后 可通过互相关在两个记录位置处得到的道并在所有边界源基础上加和而计算。如 果待反演的数据在目标区域内具有多个源和接收器（相对于在边界上具有一个或 另一个），则这是计算所需格林函数十分有效的方法。但是，对于地震数据的情况， 待反演数据的源和接收器均在目标区域中是很少的。因此，这种改进对于地震反 演问题具有非常有限的应用性。

Berkhout[7]和Zhang[8]提出反演总体上可通过反演非编码同时源而改进，该 非编码同时源在地下一些区域中被相干加和以生成一定程度上所需的波前。例如， 点源数据可随时移而加和，该时移是源位置的线性函数，从而以相对于地表的一 定特定角度生成下行平面波。这种技术可用于所有种类的反演。此方法的问题是， 源集合的相干加和必定降低数据的信息量。因此，例如，加和生成平面波去除了 地震数据中所有关于传播时间相对于源-接收器炮检距（source-receiver offset）的 信息。该信息对于更新缓慢改变的背景速度模型是至关重要的，因此Berkhout法 没有很好地进行限制。为克服这个问题，可反演数据的多个不同相干加和（例如， 多个具有不同传播方向的平面波），但是然后由于反演成本与反演的不同加和数量 成比例而失去效力。在本文中，这种相干加和的源被称为广义源。因此，广义源 可以是点源或生成具有一定所需形状的波前的点源总和。

Van Riel[9]提出通过如下进行的反演：输入地震数据的非编码叠加或部分叠加 （关于源-接收器炮检距），然后限定关于此叠加数据的成本函数，其将被优化。因 此，此公开提出利用非编码的同时源改进基于成本函数的反演。与Berkhout[6]的 同时源反演方法一样，通过此方法提出的叠加减少了待反演数据中的信息量，因 此该反演比对于原始数据得到较不充分的限制。

Mora[10]提出反演数据，该数据是宽间距的源的总和。因此，此公开提出利 用非编码的同时源模拟提高反演效力。加和宽间距的源具有这样的优势：比 Berkhout提出的相干加和保留更多信息。但是，宽间距的源的加和意味着必须用 于反演的孔隙（aperture）（反演的模型区域）必须增加以容纳所有宽间距的源。由 于计算时间与此孔隙的面积成比例，Mora法没有产生如加和的源彼此接近时所可 以实现的一样多的效力增加。

Ober[11]提出通过利用同时编码源加速地震偏移——非迭代反演的特殊情况。 在测试不同编码方法后，Ober发现，所得偏移图像具有显著降低的信噪比，这是 因为宽带编码函数必须仅大约正交的事实。因此，当加和大于16个炮点时，反演 质量不理想。由于非迭代反演开始不是非常昂贵并且由于需要高信噪比反演，因 此该技术在地球物理学工业中未被广泛实践。

Ikelle[12]提出通过同时模拟在不同时间间隔（在模拟中）激活的点源而快速 正向模拟的方法。还讨论了将这些时移的同时源模拟数据解码返回单独的模拟的 方法，该单独的模拟将已得自各个点源。然后这些解码数据可被用作任何常规反 演程序的部分。Ikelle方法的问题是，提出的解码方法将生成单独的数据，其噪声 水平与相邻源的数据之间的差异成比例。此噪声对于非横向恒定的地下模型将是 显著的，例如包含倾斜反射器（dipping reflectors）的模型。此外，此噪声将与同 时源数成比例增长。基于这些困难，如果用于反演非横向恒定的地下，Ikelle的同 时源方法可导致不可接受水平的噪声。

Krebs等在PCT专利申请公开号WO2008/042081中提出的源编码——本文将 其被赋予的全部权限引入作为参考——是反演全波场数据非常成本有效的方法。 （同时反演编码的集合的同样的方法将服务于接收器——通过源-接收器互易性或 通过编码共源数据集合中的实际接收器位置）。对于固定接收器，单个有效源仅需 进行正向和伴随计算；参见PCT专利申请公开号WO2009/117174，本文将其被赋 予的全部权限引入作为参考。如果事实是针对一般的2D采集几何记录数百炮点和 在3D考察的情况下记录数千炮点，此方法带来的计算存储（saving）相当大。在 实践中，对于最常见的场数据采集几何，固定接收器假设不严格地有效。在海洋 拖缆数据的情况下，每一个新炮点的源和接收器均发生移动。即使是在接收器位 置固定的考察中，实践通常是：不是所有接收器都“监听”每一个炮点，并且正在监 听的接收器可自炮点-至-炮点而不同。这也有悖于“固定接收器假设”。此外，由于 逻辑问题，难以记录接近源的数据，并且这意味着一般缺少炮检距附近的数据。 对于海洋和陆地考察也是这样的。这些因素均意味着，对于同时源集合，每一个 接收器位置均将漏测一些源炮点的数据。总之，在同时编码源反演中，对于给定 的同时编码集合，需要每一个炮点的所有接收器位置的数据，并且这可被称为同 时编码源反演的固定接收器假设。在WO08/042081中，当固定接收器假设不理想 时，一些公开的实施方式可比其它的好。因此，对同时编码源（和/或接收器）反 演的直接应用进行调节或调整将是有利的，这在有悖于固定接收器假设时将增强 其性能。本发明提供了如此实施的方式。Haber等[15]也描述了关于利用随机优化 方法移动同时编码源反演中的接收器的问题的方法，并将其应用于直流电阻性问 题。

发明概述

在一个实施方式中，本发明是计算机执行的测量地球物理数据的全波场反演 方法，从而确定地下区域的物理性质模型，该方法包括：（a）通过同时编码源和/ 或接收器反演利用计算机反演从测量地球物理数据选定的波至浅层时窗（time  window of arrivals），从而得到地下区域的第一物理性质模型；（b）通过迭代顺序 源反演——可仅利用测量数据的稀疏采样——利用计算机反演测量地球物理数据 或从测量地球物理数据选定的波至深层时窗，从而得到地下区域的第二物理性质 模型，其中第一物理性质模型被用作起始模型，并且源位置组被用于更新迭代顺 序源反演中的第二物理性质模型；和（c）将地下区域的第二物理性质模型输出、 显示或保存至数据存储器。

在另一实施方式中，本发明是计算机执行的测量地球物理数据的全波场反演 方法，从而确定地下区域的物理性质模型，该方法包括：（a）通过同时编码源和/ 或接收器反演利用计算机反演从测量地球物理数据选定的波至浅层时窗，从而得 到地下区域的第一物理性质模型；（b）利用第一物理性质模型以利用计算机模拟 与来自比所述浅层时窗更深处的波至相应的更长炮检距的合成数据；（c）利用计 算机反演测量地球物理数据，其中增加炮检距较长的数据，所述反演是同时编码 源和/或编码的接收器反演，从而得到地下区域的第二物理性质模型，其中所述增 加的炮检距较长的数据是较长炮检距的合成数据和在较长炮检距处的测量数据的 总和；和（d）将地下区域的第二物理性质模型输出、显示或保存至数据存储器。

附图简述

由于专利规定对颜色使用的限制，附图是彩色原图的黑白复制品。带有彩图 的本专利申请或公开复本可通过请求和支付必要的费用从美国专利商标局（U.S. Patent and Trademark Office）（美国申请号12/903,744）得到。

本发明及其优势将通过参考下文详述和附图得到更好的理解，在附图中：

图1是显示可用于同时源反演的数据窗的示意图；

图2是显示本发明方法的一个实施方式中的基本步骤的流程图，其中同时源 编码用于浅层时窗，并且稀疏的顺序源反演用于方法I中的较深窗，而在方法II 中浅层模型用于计算先被编码实施同时源反演的漏测数据道；

图3显示生成测量数据的实施例中所用的“真”速度模型；

图4显示利用图3的速度模型所计算的和利用WO2008/042081所述的二进制 编码所编码的顺序源数据的两个炮点集合，显示0-3秒的浅层时窗；

图5显示测试实施例中全波场反演的起始速度模型；

图6显示通过反演来自浅层窗（0-3秒）的低频数据（峰频率为7Hz）得自同 时源反演的模型；如图4所示，固定接收器假设对时窗0-3秒有效；

图7显示利用数据窗为0-4秒的低频数据（峰频率为7Hz）的同时源反演得 到的模型；由于固定接收器假设不再有效，在反演模型中可明显看到假象；

图8显示利用数据窗为0-5秒的低频数据（峰频率为7Hz）的同时源反演得 到的模型；由于固定接收器假设不再有效，在反演模型中可明显看到假象；

图9显示在测试实施例中通过利用数据窗为0-4秒的低频数据（峰频率为7Hz） 的顺序源反演、利用炮点间距为1.2km的稀疏源组得到的模型；由于应用顺序源， 模型与图7相比不具有任何假象；

图10显示在测试实施例中利用数据窗为0-5秒的低频数据（峰频率为7Hz） 的顺序源反演、利用炮点间距为1.2km的稀疏源组得到的模型；由于应用顺序源， 模型与图8相比不具有任何假象；

图11显示在测试实施例中利用数据窗为0-5秒的全带频数据（峰频率为40Hz） 的顺序源反演、利用炮点间距为1.2km的稀疏源组得到的模型；由于应用顺序源， 模型与图8相比不具有任何假象；

图12显示在测试实施例中利用图3的真模型计算的近炮检距和远炮检距有限 的炮点集合；最近的可用炮检距为200m，而最远的可用炮检距为5000m；

图13显示图12的近炮检距有限的炮点集合，在x=200处加载有最近可用炮 检距道；

图14示例利用图12所示的近炮检距有限的炮点集合的编码数据进行的反演， 用低频数据（峰频率为7Hz）的同时源反演得到在迭代5所示的模型；

图15示例利用图12所示的近炮检距有限的炮点集合的编码数据进行的反演， 用低频数据（峰频率为7Hz）的同时源反演得到在迭代10所示的模型；其明确显 示模型已经分叉，具有显著假象；

图16示例利用图13所示的近炮检距加载的炮点集合的编码数据进行的反演 （因为近炮检距加载是近似的，水层梯度被减弱以降低近似对同时源反演的影 响），用低频数据（峰频率为7Hz）的同时源反演得到在迭代5所示的模型；和

图17示例利用图13所示的近炮检距加载的炮点集合的编码数据进行的反演 （因为近炮检距加载是近似的，水层梯度被减弱以降低近似对同时源反演的影 响），用低频数据（峰频率为7Hz）的同时源反演得到在迭代200所示的模型。

将结合实施例实施方式对本发明进行描述。但是，在下文详述针对于本发明 的具体实施方式或具体应用方面，这仅意为示例性，并不被解释为限制本发明的 范围。相反，意为覆盖所附权利要求限定的可包括在本发明范围内的所有替代选 择、修正和等同形式。

实施例实施方式详述

本发明的一个实施方式是同时编码源反演与传统顺序源反演的混合组合。该 实施方式将同时源编码用于数据的较浅时窗并将稀疏顺序源用于反演数据的较深 部分。Krebs等[5,16]显示，编码同时源成本函数可比常规成本函数更有效地计算， 同时仍提供准确的反演。同时源成本函数在此被定义为（相对于上述方程式（2））：

$S_{\mathrm{sim}}\left(M\right)=\underset{G=1}{\overset{N_G}{\Sigma }}W\left(\delta (M,\underset{g\in G}{\Sigma }{c}_g\otimes w_g)\right)$    （方程式6）

其中，如在方程式（2）中，暗指基于接收器和时间样本的加和，并且：

定义了通过集合子组的基于集合的总和，

S_sim=同时源数据的成本函数，

G=同时广义源组，和

N_G=组数，

c_g=与各集合源特征褶积以编码集合的时间函数，这些编码函数经选择 以关于一些适当的操作大约正交，这取决于加权函数W。当W是L2-标准，适当 的操作是互相关。

方程式（6）的外加和是基于相应于集合类型（例如，共炮点集合的点源）的 同时广义源组。基于g的内加和是基于分组进行同时计算的集合。关于一些正向 建模方法，如有限差分建模，加和的广义源（基于g∈G的内总和）的正向模型计 算可与单个源的计算进行相同时间量。因此，如Krebs等[5]所示，可利用算法3极其有效地计算。

算法3–计算编码同时源成本函数的算法。

再如Krebs等[5]所示，此算法可计算S_sim(M)比从方程式（2）计算S(M)快N_g/N_G倍。

在固定接收器的几何学中，炮点集合具有完整的接收器覆盖。对于较浅窗， 如果可居于（populated）近炮检距，则其可实现。较浅窗被限定为相应于时间小 于在最长炮检距的较快模式的波至时间。编码同时源方法可用的浅层窗的尺寸取 决于近地表速度和最大可用炮检距。图1是显示可用于同时源反演的数据窗的示 意图。斜率较小的线限定较快浅层速度情况（100）的窗底部。斜率较大的线限定 较慢浅层速度情况（101）的窗底部。垂直轴是时间，水平轴是源-接收器炮检距， 选定接收器的零炮检距处于图中间，正炮检距向右，并且负炮检距向左。图1所 示示意图表明，对于慢近地表速度，可用窗（101）相对于近地表速度较快时的情 况（100）较长。

炮点集合中的负炮检距可利用互易性（201）加载。漏测的近炮检距可用估计 值——例如来自邻道（201）——加载。加载浅层窗中的近炮检距将提供完整的接 收器覆盖以应用同时源编码。另一选择是简单地排除对于一定炮点为近炮检距的 所有接收器。存在多种在近炮检距数据中加载的现有方法，如重建方法，一般用 于地表相关的多种衰减。这使得浅层窗炮点集合（202）符合同时源编码（203） 优选的固定接收器假设。因此，对于浅层窗，可生成同时源编码数据（204）和对 拟合这些数据的地下模型（206）进行反演（205）。由于计算时间明显减少几个数 量级，如必要有可能可使用更精密的模拟算法（例如，弹性全波场反演）。在用浅 层窗数据进行同时源反演后，下一步是使用浅层模型作为起始模型（207）以反演 较深窗数据（208），因为固定接收器假设对于较深窗是失效的。

本发明方法包括两种反演较深窗数据的实施方式。这些混合方法的流程图示 于图2。在图2所示的第一方法中，利用传统顺序炮点技术反演（209）较深窗数 据（208）。但是，较稀疏源组足以更新模型（210）的较深部分，并可因此有利地 使用。需要多少源将取决于模型的空间波长和最大可用的炮检距。从物理学角度 而言，地下的分辨率随深度下降。例如，在炮检距有限的数据中，较深窗的贡献 主要来自于反射而非透射。因此，取决于需要更新的模型所需的规模长度，必须 有多少顺序源的选择可与问题相关。例如，如果目标只是得到更好地解释运动学 （事件进程时间（event traveltimes））的平滑速度模型，则本发明的使用者可选择 稀疏源组用于反演。如果目标是得到高分辨率靶向模型更新，则使用者可在目标 区域上方包括更多源和在其他各处包括稀疏源。

在第二方法（图2的方法II）中，优势由如下事实形成：源自浅层窗反演的浅 层模型（206）应满足预期有时大于浅层窗的长炮检距数据（211）。因此，利用得 自同时反演的浅层模型计算较大漏测炮检距（211）的数据。此方法的优势是，在 利用正向建模得到的可用远炮检距道的情况下，固定接收器假设可适用。因此， 在可用数据连同计算远炮检距的情况下，炮点集合可被编码（212）用于较深窗 （213）反演（214），如图2所示。

本发明的混合方法不仅可用于拖缆数据，还可用于固定接收器假设失效的多 种其他采集几何。例如，海底电缆（OBC）采集一般基于插入码（patch-based）（炮 点子组的固定接收器），并且其不符合所有接收器均固定并且记录所有炮点的理想 化情况。类似于陆地采集，由于逻辑问题，固定接收器的几何难以实现。

多种类型的编码函数c_g可用于方程式（6），包括但不限于：

·线性、随机、线性调频脉冲和改良线性调频脉冲频率相关的相位编码，如 Romero等[13]所示；

·频率不相关的相位编码，如Jing等[14]所示；

·随机时移编码；

·用于远程通信的频率划分多路传输（FDMA）、时间划分多路传输（TDMA） 和编码划分多路传输（CDMA）。

取决于用途，这些编码技术中的一些将比其它的好，并且一些可组合。具体 地，利用频率相关的随机相位编码以及通过组合近源的频率不相关编码与更宽间 距源的频率相关随机相位编码已取得良好的效果。不同编码的相对优势可通过用 各组编码函数运行测试反演以确定哪种更快收敛而得到指示。

应当注意，同时编码源技术可用于多种类型的反演成本函数。具体地，其可 用于基于上述L2以外的其他标准的成本函数。其也可用于比方程式2所示成本函 数更精密的成本函数，包括正则化成本函数。最后，同时编码源方法可与任何类 型的全局或局部成本函数反演方法一起使用，包括蒙特卡罗、模拟退火、遗传算 法、进化算法、梯度线搜索、共轭梯度和牛顿法。

本发明方法也可联合多种类型的广义源技术使用，如Berkhout[7]提出的那些。 在这种情况下，不是编码不同的点源集合特征，而是编码不同合成平面波的特征。

对上述实施方式一些变型包括：

·反演每次迭代的c_g编码函数可改变。在至少一些实例中，这导致反演更 快收敛。

·在一些情况下（例如，当源采样密度大于接收器采样时），利用互易性处 理作为计算源的实际接收器和编码接收器而非源可以是有利的。

·本发明不限于单组件点接收器。例如，接收器可以是接收器阵列或其可以 是多组件接收器。

·本方法可通过优化编码而改进，以产生最高质量反演。例如编码函数可经 优化以减少成本函数中的局部最小值数。编码函数可通过对利用不同编码 函数进行的测试的人工检查或利用自动化优化程序进行优化。

·同时编码源数据的采集可导致大量地球物理数据采集成本的节约。

·对于海洋地震数据考察，从多个同时运行的海洋振动器采集编码源数据将 是十分有效的，该海洋振动器在运转中连续运行。

·如上所述，本发明的编码过程可在数据野外采集中进行，例如其中用不同 的编码函数编码多个同时运行的振动器的指示信号。在所附权利要求中， 涉及编码地球物理数据、或涉及来自编码源的地球物理数据、或涉及得到 地球物理数据的编码集合的步骤将被理解为包括得到在野外采集过程中 已经编码的数据，除非上下文明确指出编码发生在数据处理步骤。

·可以使用关于成本函数的其他限定，包括不同标准（例如，L1标准（绝 对值）替代L2标准）和其他正则化和稳固化反演的项（例如，将处罚不 平滑的模型或不稀疏的模型的项）的应用。

实施例

图3-11显示当有悖于固定接收器假设时利用本发明的混合方法反演恒定密度 声学地震数据的合成实施例。将结果与固定接收器假设有效的同时源反演进行比 较。

图3是真速度模型，即，将用于生成合成数据的速度模型。模型的水深为500 m，并且储集层的深度为3km。

图4显示两个远炮检距有限的代表性炮点集合的实施例，其用PCT专利申请 公开号WO2008/042081所述的二进制编码进行编码。由图显而易见的炮点集合的 炮检距有限的特征使得固定接收器假设无效。但是，还可看出固定接收器假设对 于0至约3秒的浅层时窗是有效的。图4显示x=800与x=2000之间的左侧炮点 的接收器不具有任何来自右侧炮点的贡献，即来自右侧炮点的能量在3秒前没有 到达那些接收器。类似地，x=2800与x=4200之间的右侧炮点的接收器在0至3 秒之间浅层时窗中不具有任何来自左侧炮点的贡献。此实施例所用的最远炮检距 是5km。数据用80m炮点间距生成，并且接收器间距为10m。本实施例利用近 炮检距数据，该近炮检距数据用模型计算得到。为生成反演的起始速度模型，利 用5-km炮检距有限的顺序炮点生成共深度点（CDP）集合，并利用Kirchoff深度 偏移展平共同的图像集合。图5所示的起始速度模型利用Kirchoff深度偏移而得到。 由于深度偏移拟合数据的运动学部分，起始模型是平滑的模型。

下一步利用同时源反演用图5的起始模型反演编码数据。由于全波场反演是 高度非线性问题，时间-频率窗一般是必需的，从而使稳定收敛的适定问题得到期 望的解决。

图6显示通过反演低频数据（峰频率为7Hz）得到的模型。由于上至3秒的 数据不违背固定接收器假设（图4显示），在同时源反演的情况下恢复的模型实现 稳定的收敛。对于上至3s的时窗数据，最大炮检距为约5km。因此反演相对于较 深部分很好地恢复（重新获得，recover）模型的较浅部分。图7和8显示在时窗 分别为4s和5s的情况下反演数据的结果。由结果明确的是，反演模型具有假象， 因为这些时窗中的编码数据不具有大的炮检距数据，因此对于模型更新计算了不 正确的梯度。由于在这些时窗中固定接收器假设被违背，所以同时源反演在模型 中生成假象。

得自浅层窗同时源反演的模型（图5）被用作较深时窗4s的顺序源反演的起 始模型。稀疏源组（用于编码数据的383个源中的20个源）用于4s窗反演。所 得模型显示在图9中。利用图9的模型作为起始模型，反演5s窗数据。对于5s 数据反演，仅使用10个顺序源。所得模型显示在图10中。对于炮检距有限的数 据，孔隙随深度而变得有限，因此稀疏源组足以得到合理的模型。利用稀疏源使 得计算作业显著减少。

图11显示用10个顺序源反演全带数据（峰频率为40Hz）的5s数据窗的结 果。

在上述实施例中，使用测量近炮检距数据。在实践中，海洋拖缆数据一般具 有漏测的近炮检距。在下一个实施例中，排除上至200m的炮检距数据。其中目 标是显示漏测近炮检距在同时反演中的影响及其被包括在编码数据中的重要性。

图12显示代表性炮点集合，其具有漏测炮检距，该漏测炮检距被编码于同时 反演。由于编码数据中存在漏测近炮检距信息，同时反演不仅在恢复的模型中生 成假象，而且还难以拟合数据。在迭代5和50恢复的的模型分别显示在图14和 15中。反演模型具有大量假象，并且该方案对于低频数据反演已经完全发散。图 15显示，由于编码数据中的漏测近炮检距信息，水层中存在大量的速度假象。为 改进该方案，可利用数据内插或正则化在近炮检距道中进行加载。这是为多种排 除问题制备数据的普遍程序。另一导致反演中假象的因素是由于直接波至导致的 漏测大能量数据。一般，直接波至的近炮检距数据具有巨大能量，因此从数据拟 合的角度而言，其对数据错配度具有大贡献。测量编码数据中的漏测信息加剧问 题并导致反演发散，造成显著假象。认识到数据拟合的这个方面，可通过屏蔽或 减弱通常与水速相关的模型浅层部分中的梯度而在一定程度上规避该问题。这在 下一个实施例中证明。

取代内插近炮检距道——其是用于多种排除的普遍程序，利用200m处的最 近可用道，并制备多份以居于（populate）漏测近炮检距道，该道显示在图13中。 然后利用二进制编码编码数据，进行同时反演。在浅层3s时窗数据的低频编码数 据的同时反演过程中，水层梯度被减弱，使得较深部分的模型更新的水层中积累 的误差不传播，并导致反演发散。

图16和17分别显示迭代5和200的反演结果，其利用在水层梯度减弱下低 频数据的同时源反演。图17清楚表明，反演是稳定的，并且恢复的模型类似于图 5的反演。如果利用数据正则化方法如radon内插内插近炮检距，则结果应进一步 改进。所示所有实施例均显示图2所描述的方法“I”。

前述申请涉及本发明的具体实施方式，其目的是为其示例。但对本领域技术 人员显而易见的是，本文所述实施方式的多种修正和变型是可能的。所有这种修 正和变型均意图在所附权利要求限定的本发明的范围内。本领域技术人员将容易 理解，在本发明的优选实施方式中，本发明方法中的至少一些步骤在计算机中进 行，即本发明是计算机实施的。在这种情况下，所得更新物理性质模型可被下载、 显示或保存于计算机存储器。

1.Tarantola,A.,"Inversion of seismic reflection data in the acoustic approximation," Geophysics49,1259-1266(1984).

2.Sirgue,L.,and Pratt G."Efficient waveform inversion and imaging:A strategy for  selecting temporal frequencies,"Geophysics69,231-248(2004).

3.Fallat,M.R.,Dosso,S.E.,“Geoacoustic inversion via local,global,and hybrid  algorithms,”Journal of the Acoustical Society of America105,3219-3230(1999).

4.Hinkley,D.and Krebs,J.,“Gradient computation for simultaneous source  inversion,”PCT Patent Application Publication No.WO2009/117174.

5.Krebs,J.R.,Anderson,J.A.,Neelamani,R.,Hinkley,D.,Jing,C.,Dickens,T., Krohn,C.,Traynin,P.,“Iterative inversion of data from simultaneous geophysical  sources,”PCT Patent Application Publication No.WO 2008/042081.

6.Van Manen,D.J.,Robertsson,J.O.A.,Curtis,A.,"Making wave by time reversal," SEG International Exposition and 75^th Annual Meeting Expanded Abstracts,1763-1766 (2005).

7.Berkhout,A.J.,"Areal shot record technology,"Journal of Seismic Exploration1, 251-264(1992).

8.Zhang,Y.,Sun,J.,Notfors,C.,Gray,S.H.,Cherris,L.,Young,J.,"Delayed-shot 3D depth migration,"Geophysics70,E21-E28(2005).

9.Van Riel,P.,and Hendrik,W.J.D.,“Method of estimating elastic and  compositional parameters from seismic and echo-acoustic data,”U.S.Patent No. 6,876,928(2005).

10.Mora,P.,"Nonlinear two-dimensional elastic inversion of multi-offset seismic  data,"Geophysics52,1211-1228(1987).

11.Ober,C.C.,Romero,L.A.,Ghiglia,D.C.,"Method of Migrating Seismic  Records,"U.S.Patent No.6,021,094(2000).

12.Ikelle,L.T.,"Multi-shooting approach to seismic modeling and acquisition,"U.S. Patent No.6,327,537(2001).

13.Romero,L.A.,Ghiglia,D.C.,Ober,C.C.,Morton,S.A.,"Phase encoding of shot  records in prestack migration,"Geophysics65,426-436(2000).

14.Jing X.,Finn,C.J.,Dickens,T.A.,Willen,D.E.,"Encoding multiple shot gathers  in prestack migration,"SEG International Exposition and70^th Annual Meeting  Expanded Abstracts,786-789(2000).

15.Haber,E.,Chung M.and Herrmann,“An effective method for parameter estimation  with PDE constraints with multiple right hand sides,”Preprint-UBC  http://www.math.ubc.ca/～haber/pubs/PdeOptStochV5.pdf(2010).

16.Jerome R.Krebs,John E.Anderson,David Hinkley,Ramesh Neelamani, Sunwoong Lee,Anatoly Baumstein,and Martin-Daniel Lacasse,“Full-wavefield  seismic inversion using encoded sources,”Geophysics74-6,WCC177-WCC188(2009).