CN101952743B - 地球物理数据反演的有效方法 - Google Patents

Abstract

一种对来自地下区域的测量的地球物理数据进行有效反演以勘探碳氢化合物的方法。测量的数据的集合(40)使用一组不等效的编码函数(30)编码(60)。之后，与单个接收器对应的每个编码的集合中的所有数据记录被求和(60)，针对每个接收器重复该步骤以产生同步的编码集合(80)。该方法使用迭代的、成本函数的局部优化来反演同步源数据的编码集合。一种伴随方法用于计算局部优化过程(100)所需的成本函数的梯度。在迭代更新之后，反演的数据产生地下区域的物理性质模型(110)，所述模型可以指示存在碳氢化合物堆积物。

Description

地球物理数据反演的有效方法

相关申请的交叉引用

本申请要求于2008年3月21日递交的美国临时申请61/070,297的优先权，并且与于2008年4月10日公开、公开号为WO2008/042081的PCT申请有关，其全文通过引用合并于此。

技术领域

本发明一般涉及地球物理勘探领域，并且更具体地涉及地球物理数据处理。具体地，本发明是一种对从多个地球物理源诸如地震源获得的数据进行反演的方法，涉及地球物理仿真或模拟，其在仿真(模拟)的一次执行中对来自很多同时活动的地球物理源的数据进行计算。

背景技术

地球物理反演[1，2]尝试找到一种地下性质模型，该模型最优解释观察到的数据并且满足地质学和地球物理的约束。存在大量公知的地球物理反演方法。这些公知的方法属于迭代反演和非迭代反演两种类别中的一种。以下内容是两种类别中的每种通常表示的定义：

非迭代反演——通过假设某种简单背景模型以及基于输入数据更新模型来实现的反演。这种方法并不将更新的模型用作另一反演步骤的输入。对于地震数据的情形，这些方法通常被称为成像、迁移、衍射层析成像或博恩(Born)反演。

迭代反演——涉及地下性质模型的反复改进的反演，从而找到一种很好地解释观察到的数据的模型。如果反演收敛，则最终的模型将更好地解释观察到的数据并且将更接近地近似实际的地下性质。迭代反演通常产生比非迭代反演更准确的模型，但是计算更为昂贵。

迭代反演通常优于非迭代反演，因为其产生更为准确的地下参数模型。遗憾的是，迭代反演在计算上如此昂贵以致于对很多感兴趣的问题应用迭代反演是不现实的。这种高计算开销是所有反演技术需要很多密集型计算仿真的结果。任何单个仿真的计算时间与要反演的源的数目成比例，并且地球物理数据中通常存在大量源。迭代反演加剧了此问题，因为必须要计算的仿真的数量与反演中迭代的次数成比例，并且需要的迭代次数通常大约是数百到数千。

地球物理中最常采用的迭代反演方法是成本函数优化。成本函数优化涉及相对模型M的成本函数S(M)的值的迭代最小化或最大化，所述成本函数S(M)是对计算的和观察到的数据之间的不吻合(或不匹配)的度量(这有时也称为目标函数)，其中计算的数据使用计算机利用当前地球物理性质模型和在给定的地球物理性质模型所表示的介质中源信号的物理学控制传播来仿真。仿真计算可以利用多种数值方法中的任意一种来实现，这些方法包括但不限于有限差分、有限元或光线追踪。模拟计算可以在频域或者时域中执行。

成本函数优化方法要么是局部的要么是全局的[3]。全局方法仅涉及计算模型群{M₁，M₂，M₃，...}的成本函数S(M)以及从模型群中选择一组一个或更多近似最小化S(M)的模型。如果需要进一步改进，则之后可以将该新选择的模型组用作产生新的模型群的基础，所述新的模型群可以相对成本函数S(M)被再次测试。对于全局方法，测试群中的每个模型可以被视为一次迭代，或者在更高水平每组测试的模型群可以被视为一次迭代。众所周知的全局反演方法包括蒙特卡洛法(Monte Carlo)、仿真退火、遗传和进化算法。

遗憾的是，全局优化方法通常收敛得很慢并且因此大多数地球物理反演是基于局部成本函数优化。算法1概述了局部成本函数优化。

1.选择初始模型，

2.关于描述模型的参数，计算成本函数S(M)的梯度，

3.搜索更新的模型，该模型是逆梯度方向中初始模型的扰动，其更好地解释了观察到的数据。

算法1——用于执行局部成本函数优化的算法。

这个过程通过将新更新的模型用作针对另一梯度搜索初始模型来进行迭代。该过程持续直到找到很好地解释所观察数据的更新模型。通常使用的局部成本函数反演方法包括梯度搜索、共轭梯度和牛顿(Newton)方法。接下来，将略为详细地解释此背景信息。

声学近似中地震数据的局部成本函数优化是普通的地球物理反演任务，并且通常说明其他类型的地球物理反演。当反演声学近似中的地震数据时，成本函数可以写成：

$S\left(M\right)=\underset{g=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{r=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{t=1}{\overset{N_t}{\mathrm{\Σ}}}W({\mathrm{\ψ}}_{\mathrm{calc}}(M,r,t,w_g)-{\mathrm{\ψ}}_{\mathrm{obs}}(r,t,w_g))---\left(1\right)$

其中：

S＝成本函数；

M＝描述地下模型的N个参数(m₁，m₂，...m_N)的向量，

g＝集合索引，

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

N_g＝集合的数量，

r＝集合内的接收器索引，

N_r＝集合中接收器的数量，

t＝迹内的时间样本索引，

N_t＝时间样本的数量，

W＝最小化标准函数(我们通常选择W(x)＝x²，其为最小平方(L2)标准)，

ψ_calc＝自模型M计算的地震压力数据，

ψ_obs＝测量的地震压力数据。

方程式1中的地震数据的集合可以是能够在地震前向建模程序的一次运行中被模拟的任意类型的集合。通常，这些集合对应于地震炮，然而这些地震炮可以比点源更为普通。对于点源，集合索引g对应于各点源的位置。对于平面波源，g将对应于不同的平面波传播方向。可以在野外获得此广义源数据ψ_obs或者可以从利用点源获得的数据中合成此广义源数据ψ_obs。另一方面，计算的数据ψ_obs通常可以在前向建模时通过利用广义源函数来直接计算。对于很多类型的前向建模，包括有限差分建模，广义源所需的计算时间大致等于点源所需的计算时间。

方程式1可以简化为：

$S\left(M\right)=\underset{g=1}{\overset{N_g}{\mathrm{\Σ}}}W\left(\mathrm{\δ}(M,w_g)\right)---\left(2\right)$

其中暗含在接收器和时间样本上的求和并且，

δ(M，w_g)＝ψ_calc(M，w_g)-ψ_obs(w_g)    (3)

通过成本函数优化反演的目标是尝试更新模型M，使得S(M)最小。这一点可以通过局部成本函数优化来实现，局部成本函数优化将给定的模型M^(k)更新为：

$M^{(k+1)}=M^{\left(k\right)}-{\mathrm{\α}}^{\left(k\right)}{\▿}_{M}S\left(M\right)---\left(4\right)$

其中k是迭代次数，α是模型更新的标度尺寸，

是不拟合函数的梯度，相对于模型参数获得。模型扰动或模型更新采用的值通过目标函数的梯度与步长α相乘计算，该值必须被重复计算。

根据方程式2，可以推导出成本函数的梯度的以下方程式：

${\▿}_{M}S\left(M\right)=\underset{g=1}{\overset{N_g}{\mathrm{\Σ}}}{\▿}_{M}W\left(\mathrm{\δ}(M,w_g)\right).---\left(5\right)$

因此，为了计算成本函数的梯度，需要单独计算每个集合对成本函数贡献的梯度，之后将这些贡献求和。因此，计算

所需的计算消耗是确定单个集合对梯度的贡献所需的计算消耗的N_g倍。对于地球物理问题，N_g通常对应于地球物理源的数量(源装置的每个位置被视为单独的源)并且大约为10,000到100,000，显著地放大了计算

的成本。

应当注意，

的计算需要计算W(δ)相对N个模型参数m_i的每个的导数。由于对于地球物理问题，N通常很大(不同参数的数量乘以模型网格单元数量，其中参数必须被赋值，单元数量通常大于一百万)，如果不得不针对每个各个模型参数执行计算，则此计算会及其耗时。

所需的是一种计算成本函数梯度的更为有效的方法，该方法不会显著降低局部成本函数优化的准确性。本发明满足了这种需求。

发明内容

在一个实施例中，本发明是计算机实现的方法，用于反演测量的地球物理数据以确定地下区域的物理性质模型，此方法包括：

(a)获得由测量的地球物理数据的两个或更多的编码集合构成的一个分组，其中，每个集合与单个的广义源相关或者利用源接收器互反性与单个接收器相关，并且其中每个集合使用选自一组不等效编码函数的不同编码函数被编码；

(b)通过对每个集合中与单个接收器对应的(或者在使用互反性的情况下与源对应的)所有数据记录求和以及针对每个不同的接收器重复进行的求和，来求和所述分组中的编码集合，产生同步的编码集合；

(c)假设所述地下区域的物理性质模型，所述模型提供遍布所述地下区域各位置处的至少一个物理性质的值；

(d)反演测量的地球物理数据，每次反演一个同步的编码集合，将假设的物理性质模型用作初始模型，并迭代更新所述模型以最小化测量模型仿真的数据与测量的地球物理数据之间的不吻合度的成本函数来产生更新的物理性质模型，其中利用成本函数相对于至少一个模型参数的梯度进行模型调整，所述梯度根据在时间上前向仿真的编码同步源数据和在时间上后向仿真的编码同步源数据的乘积的时间积分来计算。

更新的物理性质模型通常从计算机下载或保存到计算机存储装置，并且最终用来推断地下区域中存在或不存在碳氢化合物堆积物。

附图说明

参照以下详细说明和附图将更好地理解本发明及其优点，附图中：

图1是显示本发明的一个实施例的基本步骤的流程图；

图2是显示本发明的可替代实施例的基本步骤的流程图；

图3显示用于产生本文给出的示例中所使用的测量数据的基础速度模型，其中灰度色标条以单位米/秒显示速度；

图4显示来自图3中所示模型的有序源数据的一个源集合；

图5显示从图2的步骤260产生的同步编码源数据；

图6显示反演示例的初始模型(图2中的210)；

图7显示示例中有序源和编码同步源反演的L2范数模型拟合(fit)和迭代；

图8显示示例的有序源反演的第十二次迭代的结果；以及

图9显示示例中编码的同步源反演的第二十一次迭代的结果。

将结合本发明的优选实施例描述本发明。但是，就以下详细描述具体到特定实施例或发明的特定使用来说，仅仅是为了示例，而不是限制本发明的范围。相反，其意在覆盖可以包含在本发明的范围内的所有替代、修改和等价物，如所附权利要求所限定的。

具体实施方式

本发明体现了一种实现，其通过使用称为伴随方法(adjoint method)的方法计算局部优化过程所需的成本函数的梯度，之后应用该方法反演同步源数据的编码集合，从而计算编码的同步源成本函数的梯度，大大减少执行局部成本函数优化所需的计算时间。

伴随方法可以用来有效地针对所有模型参数一次执行(见接下来公式5的讨论)W(δ)相对于N个模型参数m_i中每个参数的导数的计算[1]。对于最小平方目标函数和网格化模型参数化的伴随方法可以通过以下算法概括：

1.将当前模型和集合签名w_g用作源来计算数据的前向仿真(forwardsimulation)以得到ψ_calc(M^(k)，w_g)

2.从仿真数据中减去观察到的数据，提供δ(M^(k)，w_g)，

3.将δ(M^(k)，w_g)作为源计算反向仿真(reverse simulation)(即时间向后)产生ψ_adjoint(M^(k)，w_g)

4.计算ψ_calc(M^(k)，w_g)和ψ_adjoint(M^(k)，w_g)的乘积对时间的积分以得到

算法2——使用伴随方法计算网格化模型的最小平方成本函数梯度的算法。

尽管使用伴随方法计算梯度相对于其他方法而言是有效的，但是这仍然是代价非常高的。特别是，伴随方法需要两种仿真，一种在时间上向前而一种在时间上向后，并且对于地球物理问题来说这些仿真通常是非常密集型的计算。此外，如上所述，这种伴随方法计算必须针对每个测量的数据集合单独执行，增加了计算成本N_g倍。

如参考文献[4]所述，可以通过对来自源的组合的数据进行反演而不是单独对各个源反演的方式减少所有类别的反演的计算成本。这可称为同步(或同时)源反演。已知的几种源组合类型包括：连贯(coherently)求和密集源(closely spaced source)以产生会产生某种所需形状的波前(如平面波)的有效源、求和宽距离源或在反演之前完全或部分堆叠数据。

通过反演组合的源获得的计算成本减少至少部分由以下事实抵消：组合的数据的反演通常产生较不准确的已反演模型。准确性的损失是由于信息在各个源被求和时丢失的事实，因此求和的数据不能像未求和的数据那样强地约束已反演模型。求和过程中的这个信息损失可以通过在求和之前编码每个炮记录(shot record)的方式得以最小化。组合之前的编码显著地保留同步源数据中的更多信息，并且从而更好地约束反演[4]。编码还允许密集源的组合，从而允许更多源针对给定计算区域被组合。可以和此技术一起使用各种编码方案，包括时移编码和随机相位编码。接下来简单回顾一下各种公布的地球物理同步源技术，既有编码的又有未编码的。

Van Manen[5]提出使用地震干涉测量法来加速前向仿真。地震干涉测量法通过将源置于感兴趣区域边界上的任何位置来起作用。这些源被单独建模并且波场在期望格林(Green)函数的所有位置被记录。之后，任何两个记录位置之间的格林函数能够通过交叉关联(相互关联)在两个记录位置获取的迹以及对所有边界源进行求和被计算。如果要被反演的数据具有位于感兴趣区域内的大量源和接收器(和在边界上具有一个或另一个源和接收器相反)，则这对于计算期望的格林函数是很有效的方法。但是，对于地震数据的情形，要被反演的数据的源和接收器都位于感兴趣区域内是很少见的。因此，这种改进对于地震反演问题具有很有限的适用性。

Berkhout[6]和Zhang[7]提出反演通常可以通过对未编码的同步源进行反演而得到改进，所述未编码的同步源被连贯求和以产生某地下区域内的某个期望波前(或波阵面)。例如，点源数据可以使用时移求和，时移(时间移位)为源位置的线性函数以产生在相对于平面的某特定角度的下行(down-going)平面波。这种技术可以适用于所有类别的反演。这种方法的问题是源集合的连贯求和必然减少数据内的信息量。因此，例如，为产生平面波的求和操作去除了与相对于源接收器偏移的传播时间相关的地震数据中的所有信息。此信息对于更新变化慢的背景速度模型是至关重要的，并且因此Berkhout的方法没有受到充分约束。为了克服这个问题，可以反演数据的多个不同连贯和(例如，具有不同传播方向的很多平面波)，但是之后效率由于反演成本与反演的不同和的数目成比例而损失了。我们称这种连贯求和的源为广义源。因此，广义源可以是点源或点源的和，其产生某种期望形状的波前。

Van Riel[8]提出通过对输入的地震数据进行未编码堆叠或局部堆叠(相对于源接收器偏移)的反演，之后相对将被优化的该堆叠的数据定义成本函数。因此，此公开文献提出基于利用未编码的同步源的反演改进成本函数。假设Berkhout[6]的同步源反演方法是适用的，则此方法提出的堆叠减少了将被反演的数据中的信息数量，并且由此与用原始数据约束反演相比，对反演的约束就不那么好。

Mora[9]提出了反演作为宽距离源的和的数据。因此，此公开文献建议了利用未编码的同步源仿真提高反演效率。对宽距离源求和的优点是保留的信息相对Berkhout提出的连贯和多得多。但是，对宽距离源的和意味着必须增加反演中必须使用的孔(反演的模型区域)以容纳所有宽距离源。由于计算时间与此孔的面积成比例，因此Mora的方法不会产生与求和的源彼此靠近所实现的效率提高同样的效率提高。

Ober[10]建议通过使用同步编码源加速地震迁移(一种非迭代反演的特定情形)。在测试各种编码方法之后，Ober发现得到的偏移图像已经显著地降低了信噪比，这归因于宽带编码功能必然只是近似正交的事实。因此，当对16个以上的炮求和时，反演的质量并不是令人满意的。由于开始非迭代反演的成本并非很大并且期望高信噪比反演，因此这种技术在地球物理行业没有得到广泛的实践。

Ikelle[11]建议通过同时仿真以变化的时间间隔激活(在仿真中)的点源的方式实现快速前向仿真的方法。还讨论了一种将这些时移的同步源仿真数据解码回单独的仿真的方法，单独的仿真本该从各个点源获得。之后，这些解码的数据被用作任何常规反演过程的部分。Ikelle方法的一个问题是提出的解码方法将产生分离的数据，其噪声等级与来自相邻源的数据之间的差异成比例。对于横向上非恒定的地下模型，此噪声将变得显著，例如来自包含倾斜反射体的模型的噪声。而且，此噪声将与同步源的数目成比例的增长。由于这些困难，如果将Ikelle的同步源方案用于反演横向上非恒定的地下模型，则会导致不可接受的噪声等级。

Krebs等人的文献[4]显示了编码的同步源成本函数可以比常规成本函数更有效地被计算，同时仍旧提供准确的反演。这种同步源成本函数在此被定义为(与上述方程式2相比)：

$S_{\left(\mathrm{sim}\right)}\left(M\right)=\underset{G=1}{\overset{N_G}{\mathrm{\Σ}}}W\left(\mathrm{\δ}(M,\underset{g\∈G}{\mathrm{\Σ}}{c}_g\⊗w_g)\right)---\left(6\right)$

其中，对接收器和时间样本的求和暗含在方程式(2)中，并且：

用集合的多个子分组限定对集合的求和，

S_sim＝同步源数据的成本函数，

G＝同步广义源的分组，以及

N_G＝分组数目，

c_g＝使用每个集合的源签名进行卷积以编码集合的时间函数，这些编码函数针对每个集合索引g被选为不同的(例如，随机相位函数的不同实现)，即不相等的(例如但不限于近似正交)。

方程式(6)中的外求和是对对应于集合类型的同步广义源的分组的求和(例如，共炮集合的点源)。对g的内求和是对为同步计算而分组的集合的求和。对于一些前向建模方法，诸如有限差分建模，针对已求和的广义源的前向模型的计算(对g∈G的内求和)与针对单个源的计算的执行时间量相同。因此，如Krebs等人的文献[4]中所示的，可以利用诸如算法3的方案很有效地计算

1.将作为源利用仿真器的单次运行仿真

2.使用c_g编码函数卷积每个测量的数据集合，然后对产生的编码集合求和(即

)，

3.从步骤1的结果中减去步骤2的结果。

算法3——计算编码的同步源成本函数的算法。

再次如Krebs等人的文献[4]中所示，该算法计算S_sim(M)的速度能够是方程式2中计算S(M)的速度的N_g/N_G倍。

上述算法，或者对于本技术领域的技术人员而言可想到的其变体可用于很有效地计算S_sim(M)并且因此可大大加快仅仅需要计算成本函数的反演，如全局反演方法。事实上，上述方法也可以加快很多梯度求值方法的成本函数梯度的计算。例如，梯度的有限差分计算包括模型M的成本函数的计算，并且邻近模型M+δM使其多个参数中的一个参数被扰动少量。

但是，地球物理反演典型的极大数量的模型参数暗示了如果使用常规的梯度计算方法(如有限差分)，则S_sim(M)的梯度的计算仍将是很昂贵的。本发明通过使用伴随方法计算S_sim(M)的梯度解决了此问题。因此，本发明的目标是计算其被定义如下：

${\▿}_M{S}_{\mathrm{sim}}\left(M\right)=\underset{G=1}{\overset{N_G}{\mathrm{\Σ}}}{\▿}_{M}W\left(\mathrm{\δ}(M,\underset{g\∈G}{\mathrm{\Σ}}{c}_g\⊗w_g)\right)---\left(7\right)$

在本发明中，可以利用伴随型方法(如算法4中概述的方法)非常高效地计算量

1.在仿真器的单次运行中将当前模型和同步源集合签名

用作源计算仿真源数据的前向仿真，以获得

2.使用c_g编码函数卷积每个测量的数据集合，之后对产生的编码集合求和(即

)，

3.从步骤1的结果中减去步骤2的结果，产生

4.将

用作源计算反向仿真(即时间上向后)，产生

5.计算

和

的乘积对时间的积分以获得

算法4——使用伴随型方法计算编码的同步源成本函数梯度的算法。

在实践中，本发明通常可用于根据图1所示的本发明的一个实施例的流程图执行局部成本函数反演。如在常规的迭代反演中，该过程开始于选择初始模型10，其是实际地下物理性质的近似。(在反演的后面迭代中更新的模型20替代初始模型。)在反演的下一步骤(30)中，创建编码函数：在步骤60中这些编码函数将用于编码测量的数据40以及在步骤70中编码源签名50。如下文所述并且参考文献[4]，可以在步骤30中创建很多不同类型的编码函数。步骤60和70产生编码的测量数据(针对每个接收器对源位置求和)80和编码的源签名90，这两种数据都使用步骤30中产生的相同编码函数组进行编码，对于非本领域技术人员的读者来说，参考文献[4]中对此过程作进一步的解释。在步骤100中，使用上述方程式7和算法4或者其显著变体计算梯度。在步骤110中，基于在步骤100中计算的梯度更新初始模型，通常采用向初始模型添加或减去梯度的缩小版的方式。步骤120检查反演收敛。这可通过使用方程式6和算法3为更新的模型计算编码的同步源成本函数以及将成本函数的值与来自之前迭代的值进行比较的方式执行。如果反演被判定为不收敛，则过程迭代回步骤30，在步骤110使用更新的模型20。当测试120显示迭代收敛或者另一终止条件满足时，则该过程结束(130)。

可以在方程式6和方程式7中使用很多类型的编码函数，包括但不限于：

·Romero等人的文献[12]中介绍的线性、随机、啁啾编码以及修改的啁啾频率相关相位编码；

·Jing等人的文献[13]中介绍的频率无关相位编码；

·随机时移编码；

·电信中使用的频分复用(FDMA)、时分复用(TDMA)和码分复用(CDMA)。

根据应用，这些编码技术中的一些技术将比其他技术工作得更好，并且可以组合其中的一些技术。特别是，已经使用频率相关随机相位编码获得了好的结果，并且也通过将附近源的频率无关编码和更宽距离间隔源的频率相关随机相位编码相组合获得了好的结果。可以通过使用每组编码函数运行测试反演来确定哪个收敛更快而获得不同编码的相关优点的指示。

应当注意可以针对很多类型的反演成本函数使用同步的编码源技术。尤其是，该技术可以用于基于除上述讨论的L2以外的其他范数的成本函数。该技术也可以用于比方程式2给出的一个成本函数更复杂的成本函数，包括正则化的成本函数。最后，同步的编码源方法可以与任意类型的全局或局部成本函数反演方法一起使用，所述反演方法包括蒙特卡洛法、仿真退火、遗传算法、进化算法、梯度线搜索、共轭梯度和牛顿(Newton)方法。

本发明方法还可以与各种类型的广义源技术结合使用，诸如Berkhout[6]建议的那些技术。在此情形下，可以针对不同的合成平面波编码签名，而非编码不同的点源集合签名。

上述实施例的一些变体包括：

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

·在一些情形中(例如，当源采样比接收器采样密集时)使用互反性将实际的接收器作为计算源对待并且编码接收器而不编码源可能是有益的。

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

·可以通过优化编码以产生最高质量反演来改进本发明方法。例如，编码函数可以被优化以减少成本函数中局部最小值的数目。编码函数可以通过使用不同的编码函数执行的测试的手动检查或者使用自动的优化程序被优化。

·同步编码源数据的获取可导致显著的地球物理数据获取成本节省。

·对于海洋地震数据勘查，自同时运转的海上振动器获得编码的源数据将是很高效的，所述海上振动器在运动的同时连续运转。

·可以使用成本函数的其他定义，包括不同范数的使用(例如，L1范数(绝对值)而非L2范数)以及正则化且稳定化反演的额外项的使用(例如，将惩罚或补偿(penalize)不平滑的模型或者不稀疏的模型的项)。

图2是本发明方法用于局部成本函数反演的可替代实施例的流程图。该流程图基本类似于图1中的流程图(相应框的标记数字加上了200)，除了以下变化：

·在步骤230中，编码函数可以在迭代之间被改变。

·添加了步骤305。在步骤305中，来自步骤300的梯度使用任何已知的方法被调节，并且为了获得更好的收敛，调节的梯度用于计算共轭梯度。在(局部)优化方案中，算法的收敛可以通过选择下降方向(梯度)作为过去的和目前的下降方向的线性组合而得到改进，这是众所周知的。这被称为共轭梯度方法。调节(或预调节)是一种改进收敛的经典公知策略。

·在步骤310中，在共轭梯度方向的线搜索用于更新模型。线搜索是标准优化策略的部分并且由以下步骤构成：在当前点开始、选择求值成本函数值的方向、沿着选择的方向移动以及在每个新位置重新对成本函数求值，以及选择提供成本函数最小值的点。该方向(在此情形下)是梯度或共轭梯度。

此外，在本发明方法的此特定说明性实施例中，在反演过程中作出以下选择：

1.测量数据的输入地震集合240是共同的点源集合。

2.编码签名230被选择为来自Romero等人的文献的随机相位签名。此签名可以通过仅生成由多个时间样本构成的序列的方式获得，所述时间样本是统一的伪随机序列。

3.在步骤260和270中，集合和签名通过将集合中的每个迹与集合的编码签名卷积的方式被编码。

4.在步骤300和320中，使用时空域中的有限差分地震建模码来执行仿真。

5.在步骤300和320中，使用L2范数计算成本函数。

示例

图3至9呈现了使用本发明反演恒定密度声地震数据的合成示例，并和常规的有序源方法进行对比。对于此示例，使用图2中描述的实施例，但是去除了步骤305。包含步骤305相对于此示例改进结果，因此，这是本发明优点的稳妥示例。

图3是用于为此示例产生测量的数据(图4)的基础速度模型。根据此模型仿真的测量数据的成功反演应当看起来像图3的限带版。

图4显示了根据图3中的模型仿真的测量数据的一个源集合。共128个源使用位于模型表面附近的源被仿真(其中图4是第一个)并且在该模型内均匀地展开。1024个地震接收器也在该模型内均匀地展开。同组接收器用于每个源位置(即接收器不随源移动)。

图5示出了产生于图2的步骤260的同步编码源数据。编码函数是3第二长随机相位函数。所有128个源使用不同的随机相位函数进行编码，之后进行求和以产生此同步编码源。请注意，图5中编码的同步源数据与图4比较看起来像随机噪声。

图6示出了针对此反演示例用作初始模型(图2的210)的模型。该模型通过使图3中显示的基础模型过度平滑的方式产生。平滑被选择为足够大以去除所有反射和小比例的速度异常，同时保留模型的平滑背景。

之后，图4和图5中的数据均被反演。当然，本发明用于反演图5中的同步编码源数据。应当注意，由于有序源数据中存在128个源集合，因此这些数据的反演每次迭代必须花费的计算耗费是反演编码的同步源数据所需的计算耗费的约128倍。128倍归因于以下事实：必须针对每个有序源单独计算梯度，并且之后这些梯度被求和以产生总的有序源梯度。另一方面，仅仅自一个集合计算编码的同步源数据的梯度。

图7示出了两种反演的模型拟合和迭代次数。有序源反演(黑色方块)以大约两倍于编码的同步源反演(白色方块)的速度收敛。(编码的同步源反演的迭代21的模型拟合约等于有序源反演的迭代12的模型拟合。)然而，由于同步源反演的每次迭代仅使用计算消耗量的1/128，因此编码的同步源反演的效率是有序源反演的64倍。

图8和图9示出了两种反演结果。第12次有序源反演和第21次编码的同步源反演被选择用于显示，因为其具有几乎相同的模型拟合。请注意，这些反演看起来都类似于图3中基础模型的限带版。还要注意除了添加少量噪声之外，编码的同步源反演看起来很像有序源反演。此噪声通过对反演结果进行图像处理或者通过向反演施加调节(图2的步骤305)而易于去除。

上述应用的目的是说明本发明的特定实施例。但是，对于本领域技术人员来说，显而易见的是对本文描述的实施例的很多修改和变化是可能的。所有这种修改和变化都处于本发明的范围之内，如所附权利要求限定的。本领域技术人员将易于认识到在本发明的优选实施例中，本发明方法的至少一些步骤在计算机上执行，即本发明是计算机实现的。在这种情形下，得到的更新的物理性质模型可以被下载或保存到计算机存储装置。

参考文献：

1.Tarantola，A.，″Inversion of seismic reflection data in the acousticapproximation，″Geophysics 49，1259-1266(1984).

2.Sirgue，L.，and Pratt G.″Efficient waveform inversion and imaging：A strategyfor selecting temporal frequencies，″Geophysics 69，231-248(2004).

3.Fallat，M.R.，Dosso，S.E.，″Geoacoustic inversion via local，global，and hybridalgorithms，″Journal oj′the Acoustical Society oj′America 105，3219-3230(1999).

4.Krebs，J.R.，Anderson，J.A.，Neelamani，R.，Hinkley，D.，Jing，C，Dickens，T.，Krohn，C，Traynin，P.，″Iterative inversion of data from simultaneous geophysicalsources，″PCT Publication No.WO2008/042081.

5.Van Manen，D.J.，Robertsson，J.O.A.，Curtis，A.，″Making wave by timereversal，″SEG International Exposition and 75th Annual Meeting ExpandedAbstracts，1763-1766(2005).

6.Berkhout，A.J.，″Areal shot record technology，″Journal of Seismic Exploration 1，251-264(1992).

7.Zhang，Y.，Sun，J.，Notfors，C，Gray，S.H.，Cherris，L.，Young，J.，″Delayed-shot3D depth migration，″Geophysics 70，E21-E28(2005).

8.Van Riel，P.，and Hendrik，W.J.D.，″Method of estimating elastic andcompositional parameters from seismic and echo-acoustic data，″U.S.Patent No.6,876,928(2005).9.Mora，P.，″Nonlinear two-dimensional elastic inversion ofmulti-offset seismic data，″Geophysics 52，1211-1228(1987).

10.Ober，C.C，Romero，L.A.，Ghiglia，D.C，″Method of Migrating SeismicRecords，″U.S.Patent No.6,021,094(2000).

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

12.Romero，L.A.，Ghiglia，D.C，Ober，C.C，Morton，S.A.，″Phase encoding ofshot records in prestack migration，″Geophysics 65，426-436(2000).

13.Jing X.，Finn，C.J.，Dickens，T.A.，Willen，D.E.，″Encoding multiple shotgathers in prestack migration，″SEG International Exposition and 70^th AnnualMeeting Expanded Abstracts，786-789(2000).

Claims (19)

1.一种用于反演测量的地球物理数据以确定地下区域的物理性质模型的计算机实现的方法，所述方法包括：

(a)获得由所述测量的地球物理数据的两个或更多编码集合构成的一个分组，其中，每个集合与单个的广义源相关或者利用源接收器互反性与单个接收器相关，并且其中每个集合使用选自一组不等效的编码函数的不同编码函数被编码；

(b)通过对每个集合中与单个接收器对应的或者在使用互反性的情况下与源对应的所有数据记录求和以及针对每个不同的接收器重复该求和，来求和所述分组中的编码集合，产生同步的编码集合；

(c)假设所述地下区域的物理性质模型，所述模型提供遍布所述地下区域各位置处的至少一个物理性质的值；

(d)反演所述测量的地球物理数据，每次反演一个同步的编码集合，将假设的物理性质模型用作初始模型，并迭代更新所述模型以将测量仿真模型的数据与所述测量的地球物理数据之间的不吻合度的成本函数最小化来产生更新的物理性质模型，其中利用所述成本函数相对于至少一个模型参数的梯度进行模型调整，所述梯度根据在时间上向前仿真的编码的同步源数据和在时间上向后仿真的编码的同步源数据的乘积的时间积分来计算；以及

(e)下载所述更新的物理性质模型或将其保存到计算机存储装置。

2.根据权利要求1所述的方法，其中反演测量的地球物理数据的同步编码集合包括

(i)用计算机仿真与所述测量的数据的同步编码集合对应的同步编码集合，利用所述假设的物理性质模型，其中所述仿真使用源签名，所述源签名被编码所采用的编码函数与编码所述测量的数据的同步编码集合使用的编码函数相同，其中整个同步编码集合在单次仿真操作中被仿真；以及

(ii)计算成本函数，其测量所述测量的数据的同步编码集合与所述仿真的同步编码集合之间的不吻合度。

3.根据权利要求1所述的方法，进一步包括将不同的不等效的编码函数组用于所述模型的所述迭代更新的至少一次迭代。

4.根据权利要求1所述的方法，其中反演所述测量的地球物理数据包括：

(i)在单次仿真运行期间，使用当前物理性质模型M并将同步源集合签名用作源计算所述编码的同步源数据的前向仿真以获得ψ_calc，其中编码所述同步源集合签名所使用的编码函数c_g与编码所述测量的数据使用的编码函数相同；

(ii)通过从步骤(i)的结果中减去所述测量的地球物理数据的同步编码集合来计算所述成本函数δ；

(iii)将δ用作仿真源，计算反向仿真，即，时间上向后，产生ψ_adjoint；

(iv)计算ψ_calc和ψ_adjoint的乘积对时间的积分以获得所述成本函数的梯度；以及

(v)将所述成本函数的所述梯度用于调整和更新模型M。

5.根据权利要求1所述的方法，其中所述测量的数据的编码集合通过将来自所述集合的所有迹与对应的编码源签名在时间上卷积而被编码，所述编码源签名是所述集合的源函数与为所述集合选择的所述编码函数的卷积。

6.根据权利要求1所述的方法，其中所述测量的数据的两个或更多编码集合通过从地球物理勘查获得数据的集合而获得，其中数据从多个同步运行、唯一编码的源装置获取。

7.根据权利要求1所述的方法，其中所述测量的地球物理数据来自地震勘查。

8.根据权利要求7所述的方法，其中广义的地震源或者是所有点源或者是所有平面波源。

9.根据权利要求2所述的方法，其中所述测量的地球物理数据包括每个源激活的测量的或算出的签名并且在仿真操作中使用的编码的源签名是通过将所述测量的或算出的源签名与编码函数在时间上卷积获得的签名，所述编码函数和用于步骤(a)中编码所述对应的测量集合的编码函数相同。

10.根据权利要求1所述的方法，其中所述编码函数的类型选自由以下项构成的组：线性、随机相位、啁啾、修改的啁啾、随机时移编码以及频率无关相位编码。

11.根据权利要求1所述的方法，其中所述编码函数针对一些源是一种类型而针对其他源是另一类型，每种类型选自由以下项构成的组：线性、随机相位、啁啾、修改的啁啾、随机时移编码以及频率无关相位编码。

12.根据权利要求1所述的方法，其中所述编码函数相对使用的所述成本函数被优化。

13.根据权利要求1所述的方法，其中步骤(d)中的所述时间上向前或向后仿真操作使用有限差分、有限元或有限体积仿真码执行。

14.根据权利要求7所述的方法，其中所述物理性质模型是地震波速率、地震弹性参数、地震各向异性参数或地震非弹性参数的模型。

15.根据权利要求1所述的方法，其中局部成本函数优化方法被用于更新所述模型，所述局部成本函数优化方法包括梯度线搜索、共轭梯度或牛顿方法。

16.根据权利要求1所述的方法，其中所述成本函数是L1范数成本函数或L2范数成本函数并且所述成本函数包括正则化项。

17.根据权利要求1所述的方法，其中所述不等效的编码函数是基本正交函数。

18.一种用于反演来自地下区域的测量的地球物理数据以勘探碳氢化合物的方法，所述方法包括：

(a)获得由所述测量的地球物理数据的两个或更多编码的集合构成的一个分组，其中，每个集合与单个的广义源相关或者利用源接收器互反性与单个接收器相关，并且其中每个集合使用选自一组不等效的编码函数的不同编码函数被编码；

(b)通过对每个集合中与单个接收器对应的或者在使用互反性的情况下与源对应的所有数据记录求和以及针对每个不同的接收器重复该求和，来对所述分组中的编码集合求和，产生同步的编码集合；

(c)假设所述地下区域的物理性质模型，所述模型提供遍布所述地下区域的各位置处的至少一个物理性质的值；

(d)反演所述测量的地球物理数据，每次反演一个同步的编码集合将假设的物理性质模型用作初始模型，并迭代更新所述模型以最小化测量模型仿真的数据与测量的地球物理数据之间的不吻合度的成本函数来产生更新的物理性质模型，其中利用所述成本函数相对于至少一个模型参数的梯度进行模型调整，所述梯度根据在时间上向前仿真的编码的同步源数据和和在时间上向后仿真的编码的同步源数据的乘积的时间积分来计算；以及

(e)使用所述更新的物理性质模型来勘探所述地下区域中的碳氢化合物。

19.一种从地下区域开采碳氢化合物的方法，其包括：

(a)执行所述地下区域的地球物理勘查；

(b)获得物理性质模型，所述模型已使用权利要求1所述的方法构建，权利要求1所述的方法通过引用合并于此；

(c)使用所述物理性质模型识别所述地下区域中的含碳氢化合物区；

(d)向所述区内钻井以及从所述井开采碳氢化合物。

Images