使用频谱整形的全波场反演的收敛速度的制作方法
【专利摘要】本发明涉及一种加速地震数据（106）的迭代反演以获取地下模型（102）的方法，其使用局部成本函数最优化。在每次迭代中已更新模型的频谱被控制以匹配地下区域的已知或估计的频谱，优选为地下P阻抗的平均振幅频谱。这种控制通过对源微波（303）和数据（302）应用频谱整形滤波，或通过对成本函数的梯度（403）应用可随时间变化的滤波完成。源微波的振幅频谱（在滤波之前）应满足D(f)=fIP(f)W(f),其中f是频率，D(f)是地震数据的平均振幅频谱，且IP(f)是地下区域（306，402）中P阻抗的平均振幅频谱或其近似。
【专利说明】使用频谱整形的全波场反演的收敛速度
[0001]相关领域的交叉参考
[0002]本申请要求2011年3月30日提交的美国临时专利申请N0.61/469，478，名为Improving Convergence Rate of Full Wavefield Inversion Using Spectral Shaping,和2011年7月15日提交的美国临时专利申请N0.61/508，440，名为Convergence Rate ofFull Wavefield Inversion Using Spectral Shaping的优先权,它们全部包含在此以供参考。
【技术领域】
[0003]本发明通常涉及地球物理勘探领域且更具体地涉及地震数据的反演，这是个广义的术语用于指基于所记录的地震数据简历地下模型的处理。特别地，本发明是使用频谱整形提升全波场反演收敛速度的方法。术语全波场反演(“FWI”)用于指一种类型的反演方法，其旨在产生以确切的量化形式全面解释所记录的地震数据的地下模型:基于作为反演结果的地下，合成地震数据的精确仿真应该高度匹配真实地震数据。
【背景技术】
[0004]地球物理反演意图找到最优解释所观察的数据且满足地理和地球物理约束的地下性质模型。存在大量已知的地球物理反演方法。这些为人熟知的方法可划入两类之一，迭代反演和非迭代反演。下面是该两类的每类通常所表达的意思的定义:
[0005]非迭代反演——通过假设一些简单的背景模型且基于输入数据更新该模型完成的反演。该方法不使用更新过的模型作为对反演的另一步的输入。对于地震数据的例子来说，这些方法通常是指成像、偏移、衍射层析成像或博恩反演。
[0006]迭代反演一涉及对地下性质模型的重复改善的反演，以便所发现的模型能很令人满意地解释所观察到的数据。如果`该反演收敛，那么最终的模型将更好解释所观察到的数据且将更加高度近似真实的地下性质。迭代反演通常产生的模型比非迭代反演产生的更加精确，但是计算过程也要昂贵得多。
[0007]在地球物理学中最常用的迭代反演方法是成本函数最优化。对于模型M，成本函数最优化涉及成本函数S (M)值的迭代最小化或最大化，该函数值是对计算数据和观察数据不匹配的度量(有时其也指目标函数)，其中在计算机上使用当前的地球物理性质模型仿真计算数据而介质中震源信号的物理组织传播(physics governing propagation)是用给定的地球物理性质模型表示的。可通过一些数字方法的任何一种方法完成仿真计算，但该仿真计算不限于有限的差、有限的元素或射线轨迹。该在频率或时间域上执行该仿真计算。
[0008]成本函数最优化可是局部的，也可是总体的。总体法简要涉及计算模型群(apopulation of models) (M17M27M3,…}的成本函数S(M)以及从群中选择近似最小化S(M)的一个或更多模型的集合。如果想要进一步的改善，然后可以该新选择的模型集合为基础产生可再次经测试与成本函数S(M)相关的新模型群。对于总体法，测试群中的每个模型可认为是迭代，或在更高水平上，每个所测试的群的集合可认识是迭代。为人熟知的总体反演方法包括蒙特卡罗，仿真退火，遗传和进化算法(Monte Carlo, simulated annealing,genetic and evolution algorithms)。
[0009]不幸的是，总体最优化方法通常收敛极慢且因此大多数地球物理反演是基于局部成本函数最优化。算法1总结了局部成本函数最优化。
1.选择起动模型
2.与描述所述模型的参数相关的成本函数S(M)的梯度
[°01°] 3,在更好解释观察数据的负梯度方向上搜索作为自动模型的扰动的更新模型
[0011]算法1—执行局部成本函数最优化的算法
[0012]通过使用新的更新模型作为另一个梯度搜索的起动模型迭代该过程。这种处理继续直到所找到的更新模型令人满意地解释观察数据。通常使用的成本函数反演方法包括梯度搜索、共轭梯度和牛顿方法。
[0013]非常常见的成本函数是真实和仿真地震道(seismic trace)差平方的和。对于这种情形而言，通过2个波场的交叉相关计算该梯度，如图1中典型全波场反演工作流程所示。从源微波的估计开始(101 )，接着初始地下模型(102)，我们通过从震源到检波点位置向前传播波(104)产生仿真地震数据(103)。通过从真实地震数据(106)减去(110)仿真数据形成数据残留(105)。这些残留然后被向后传播给地下模型(107)且与震源波场交叉相关，该波场是通过从震源位置到地下的向前传播(108)产生的。该交叉相关的结果是梯度(109)，基于该梯度地下模型被更新。在新的更新模型上重复这种处理，直到仿真和真实地震数据之间的差是可接受的。
[0014]对于不同的成本函数而言，梯度的计算可是不同的。图1中工作流程的基本元素仍是相当常见的。有关本发明的关键思想可根据使用替代成本函数和梯度计算的情形做出琐碎的改变。
[0015]相比非迭代反演通常优选迭代反演，这是因为迭代反演生成更精确的地下参数模型。不幸的是，迭代反演的计算是如此的昂贵以至于将其应用到所感兴趣的许多问题上时不切实际的。这种高计算费用是所有反演技术要求许多计算密集型仿真的结果。任何单次仿真的计算时间与要反演的源的数量成比例，且通常在地球物理数据中存在大量的源，其中在前面所用的术语“源”指源装置的致动位置(activation location of a sourceapparatus)0在迭代反演中该问题会恶化，因为需要计算的仿真的数量与反演中迭代的数量成比例，且所要求的迭代数量通常约有几百到几千个。
[0016]减少全波场反演的计算成本是使所述方法对于场规模3D应用实用的关键要求，尤其当要求高分辨率时(例如对于储层表征而言)。许多建议的方法依赖同时仿真源的思想，要么编码(例如 Krebs 等人，2009 ;Ben_Hadj-Ali 等人，2009 ；Moghaddam 和 Herrmann，2010)要么连贯加总(coherently summed)(例如 Berkhout，1992 ;Zhang 等人，2005 ；VanRiel and Hendrik，2005)。基于编码同时仿真(encoded simultaneous simulation)的反演方法常常受到污染反演结果的串扰噪声的影响且通常被数据釆集配置(对于几种方法而言，一个要求是在静止的检波点上记录数据)所限制。基于连贯加总的方法通常导致信息的损失。虽然如此，这两种类型的方法都是有用的且是正在进行的研究的对象。
[0017]一种不同的减少全波场反演的计算成本的方法是通过减少要求收敛的迭代的数量，且这是本发明的目的。本发明不受上面提到的方法的典型限制的影响，但是它不排除他们的用法。事实上，原则上它可与上面提到的同时源方法(the simultaneous-sourcemethods)的任何一种联合使用，从而潜在增加节约下的计算成本。
【发明内容】
[0018]在一个实施例中，本发明是计算机实施的方法，其用于加速地震数据迭代反演的收敛以获取地下区域中一个或更多物理参数的模型，该方法包括使用局部成本函数最优化，其中更新假设的或当前的模型以减少地震数据和模型仿真数据之间的不匹配，其中在第一迭代中控制所更新的模型的频谱且然后以匹配关于地下区域的已知或估计的频谱。
【专利附图】
【附图说明】
[0019]参考下面详细的描述和附图将更好地理解本发明和它的优势，其中:
[0020]图1是示出全波场反演基本步骤的流程图；
[0021]图2是将本发明的整形滤波应用到源微波和地震数据上的影响的示意描述；
[0022]图3是示出本发明方法的实施例的基本步骤的流程图，该方法涉及将频谱整形滤波应用到输入数据和源微波上。
[0023]图4是示出本发明方法的实施例的基本步骤的流程图，该方法涉及将频谱整形滤波应用到成本函数的梯度上。
[0024]图5是示出图4的实施例的基本步骤的流程图，图4的实施例拓展到多参数反演的情形。
[0025]图6到8示出在10次迭代(图6)、50次迭代(图7)和100次迭代(图8)后全波场反演(“FWI”)的收敛。
[0026]图9到10示出I次迭代(图9)和4次迭代(图10)后FWI迭代的收敛，使用图6到8中的实例，将本发明的整形滤波应用到地震数据和源微波上。
[0027]图11A-11D示出在用本发明方法进行频谱整形之前和之后测量的炮点道集和仿真的炮点道集的合成实例。
[0028]图12示出对来自图1IA-1ID的频谱整形的交叉相关成本函数的影响。
[0029]将联系实例实施例描述本发明。然而，就下面详细的描述是针对本发明的特定实施例或特定用途来说，这倾向于仅是描述性的，且不解读成限制本发明的范围。相反，本发明倾向于覆盖所有可包括在本发明范围内的替代、变形或等价物，如附属权利要求所定义的。本领域的技术人员将很快认识到在本发明方法的实际应用中，该方法须按照这里所教授的，经编程在计算机上实施。
【具体实施方式】
[0030]本发明方法的关键思想是基于这样的假设:对地下频谱的合理估计逻辑上是已知的。如果是这种情形，可通过保证从最初的迭代开始反演产生具有想要的频谱的地下模型来大幅减少收敛所要求的迭代次数。直观可见，这暗示不需要将计算花费在大多改变地下模型的频谱的迭代上，且结果是反演以更快的速度收敛到最终答案。
[0031]由于该思想是有意义的且实际的，须解答下面的问题:
[0032](I)通常可假设可获得对地下频谱的好的估计吗？
[0033](2)怎样才能保证反演结果具有想要的频谱且，特别地，怎样从最初的迭代中获得这种结果？
[0034]在下面两个部分提供这两个问题的答案。
[0035]估计地下模型的频谱
[0036]在关于非迭代反演的论文中，Lancaster和 Whitcombe (2000)、Lazaratos (2006)、Lazaratos和David (2008)以及Lazaratos和David (2009)引入这种思想，那就是有非迭代反演产生的模型应具有平均来说与地球地下的频谱相似的频谱，如井记录(well log)测量的一样。(这里可交替使用术语振幅频谱和频谱以值振幅对频率。)对于任意给定的地区，可通过平均本地井中记录的记录曲线的频谱推导出该目标频谱。理论上，用于正入射的反射数据的合适记录曲线是P阻抗。现实中，已经观察到大多数记录曲线的平均频谱是高度相似的。事实上，各种地理位置、深度和沉积环境的典型井记录频谱是高度相似的，所以 反演目标频谱的常规形式是健壮的且很好定义。由于井记录频谱的稳定性和健壮性，上述的公开所概述的概念可广泛用于非迭代地震反演，甚至当不能进行本地井控制。
[0037]控制反演结果的频谱
[0038]通常由几个参数表征可使用地震反演估计的地下(例如P阻抗、S阻抗、密度、P速率等)。这些不同参数的频谱基本上可不同。因此，当一个涉及地下模型的频谱时，所指的参数需要进行特指。首先处理单参数反演的情形，其中仅根据压缩(P)阻抗(阻抗是速率和密度的乘积)描述该地下模型。这是反演的一般应用，因为地下反射系数大多取决于P阻抗的变化。然后我们描述该方法是怎么拓展到多参数反演。
[0039]反演结果的频谱与地震数据的频谱以及源微波的频谱相关。尤其对于地震反射数据而言，该卷积模型表明给定地下的地震反射反应可通过地震微波的卷积和地球的反射系数来计算。假设弱散射，可示出对于正入射，反射系数函数可简单计算成P阻抗的导函数。对于斜入射而言，反射系数的计算涉及其他的弹性参数，但是这不会根本上改变这里提出的方法的概念。在频率域中，描述卷积模型的基本公式是:
[0040]D(f) =flp(f)ff(f)(I)
[0041]其中f是频率，D(f)是地震数据的平均振幅频谱，ff(f)是地震微波的振幅频谱而Ip(f)是地下P阻抗的平均振幅频谱。在时间域中计算P阻抗的导数相当于在频率域中乘以i2Jif。在当前讨论中，为了简便省略了 2 π因子，因为它不影响计算或方法的实施。因子i也被省略，因为我们仅处理振幅频谱。
[0042]刚刚所讨论的内容暗示，为了使最终的反演结果具有频谱Ip (f)，有必要使用微波，其频谱w(f)通过等式(I)和数据频谱D (f)相关联。尽管理论上等式(I)仅在等式中的振幅频谱是关于特定参数—P参数时才有效，从经验中已经观察到不同弹性参数的频谱通常是非常相似的。下面所描述的是反演问题怎样再次用公式表达以便反演从最初的反演中产生具有频谱Ip(f)的模型。
[0043]FWI中的模型更新通常计算成缩放版本的目标函数的梯度，其与模型参数有关。对于常见的L2 (最小平方)目标函数，可通过交叉相关向前传播的震源波场和向后传播的残留波场来计算梯度。向前传播的震源波场的频谱与输入微波的频谱w(f)成比例。对于第一次迭代，给定典型的反演起动模型是非常光滑的且不会产生反射，数据残留基本上等于记录数据，且因此向后传播的数据残留的频谱与D(f)成比例。因此，梯度的频谱G(f)等于2个交叉相关波场的频谱(W(f)和D(f))的乘积，进一步乘以依赖频率的因子A(f)，该因子取决于要解决反演问题的特性(例如2D对3D、声学对弹性反演、更新中的弹性参数等)。要么从理论上(例如对于2D恒密度声学反演而言，A (f) = (f1/2)，要么从实验上通过计算梯度的频谱并将它和已知的频谱W(f)和D(f)的乘积进行比较推导出所述因子。所以我们可写出:
[0044]G(f)=A(f)ff(f)D(f)(2)
[0045]让我们假设地球阻抗的频谱Ip(f)逻辑上是已知的且我们有G(f)=Ip(f)。通常这不是正确的。我们仍然可通过应用整形滤波H(f)到输入微波和数据上适本地将原始的反演问题变成新问题。新整形微波Ws(f)和整形数据Ds(f)与原始微波和数据频谱相关，关系如下:
[0046]Ws (f) =H(f)ff(f)
[0047]Ds (f) =H(f)D(f) (3)
[0048]对使用微波W(f)匹配原始数据D(f)的模型进行反演等价于对使用整形微波Ws(f)匹配整形数据Ds(f)的模型进行反演。
[0049]类似于等式(2)，我们现在写出整形梯度Gs (f):
[0050]
【权利要求】
1.一种用于加速地震数据迭代反演的收敛以获取地下区域内一个或更多物理参数的模型的计算机实施方法，其包括使用局部成本函数最优化，其中假设或当前模型经更新以减少所述地震数据和模型仿真数据之间的失配，其中在第一迭代和随后的迭代中控制已更新模型的频谱以为地下区域匹配已知或估计频谱。
2.根据权利要求1所述的方法，其中所述一个或更多物理参数是P阻抗，且其中在每次迭代中所述已更新模型的频谱经控制通过在第一迭代中对用于产生所述模型仿真数据的所述地震数据和源微波应用频谱整形滤波器以匹配所述已知或估计频谱，其中在应用所述频谱整形滤波器之前选择所述源微波以具有满足D (f) =flp (f) W (f)的频谱W (f)，其中D (f)是所述地震数据的平均频谱，f是频率，而Ip(f)是地下区域内P阻抗的平均频谱或通过地下区域不同弹性参数的频谱近似，且其中所述地下区域的已知或估计频谱是Ip(f)。
3.根据权利要求2所述的方法，其中所述频谱整形滤波器满足 H(f) = {f1/2IP(f)}/{A1/2(f)D(f)} 其中A(f)是依赖频率的因子，经确定以便 G(f)=A(f)ff(f)D(f) 其中G(f)是与所述一个或更多物理参数的一个参数有关的所述成本函数的梯度的频-1'TfeP曰。
4.根据权利要求1所述的方法，其中更新所述假设或当前模型包括计算与所述一个或更多物理参数的一个参数有关的所述成本函数的梯度，且然后至少在第一迭代中对所述梯度应用频谱整形滤波器，且 其 中使用具有满足D(f)=fIp(f)W(f)的频谱W(f)的源微波产生模型仿真数据，其中D (f)是地震数据的平均频谱，f是频率，且Ip (f)是地下区域内P阻抗的平均频谱或通过地下区域不同弹性参数的频谱近似，且其中所述已知或估计频谱是Ip(f)。
5.根据权利要求4所述的方法，其中在整形G(f)之前，所述频谱整形滤波器H(f)是通过将所述一个或更多物理参数的一个参数的所述地下区域内已知或估计频谱S (f)除以所述梯度的频谱确定的。
6.根据权利要求4所述的方法，其中所述频谱整形滤波器随时间变化。
7.根据权利要求1所述的方法，其中所述地下区域的已知或估计频谱通过平均来自地下区域的井记录的频谱获得。
8.根据权利要求7所述的方法，其中同时获取至少两个物理参数的模型，并且每个物理参数的已知或估计频谱是通过平均测量或对应于该物理参数的井记录的频谱获得。
9.根据权利要求1所述的方法，其中同时获取至少两个物理参数的模型，且在所述第一迭代及随后的迭代中每个已更新模型的频谱经控制而匹配所述地下区域内相应物理参数的已知或估计频谱。
10.根据权利要求9所述的方法，其中更新特定物理参数(用指数i表示)的假设或当前模型包括计算与所述成本函数的特定物理参数有关的梯度(％)，并然后对所述梯度应用频谱整形滤波器，其中所述频谱整形滤波器被修改以适用于在所述地下区域内的特定物理参数和其已知或估计频谱Si (f)。
11.根据权利要求10所述的方法，其中在整形前通过将Si(f)除以所述梯度的频谱确定所述频谱整形滤波器。
12.根据权利要求10所述的方法，其中所述频谱整形滤波器随时间变化。
13.根据权利要求2所述的方法，其中根据下面的标准推导所述频谱整形滤波器:在对用于产生模型仿真数据的所述地震数据和源微波应用频谱整形滤波器后，与所述一个或更多物理参数的一个参数相关的所述成本函数的梯度频谱应匹配所述地下区域的已知或估计频谱。
14.根据权利要求1所述的方法，其中从由P阻抗、S阻抗、密度、P速率和S速率组成的组中选择所述一个或更多物理参数。
15.根据权利要求1所述的方法，其中使用具有频谱W(f)的地震源微波产生所述模型仿真数据，在比例常量范围内W(f)满足下面等式:D(f)=flp(f)ff(f)其中f是频率，D(f)是所述地震数据的平均频谱，且Ip(f)是所述地下区域内P阻抗的平均频谱或基于另一个弹性参数的频谱的近似。
16.根据权利要求15所述的方法，其中D(f)和Ip(f)是平均频谱，意味着在所述地下区域内被平均。
17.一种计算机程序产品，其包括具有植入计算机可读程序代码的永久计算机可用介质，所述计算机可读程序代码适于被执行以实施用于反演地震数据从而获取地下区域内的一个或更多物理参数的模型的方法，所述方法包括使用具有局部成本函数最优化的迭代反演，其中假设或当前的模型被更新以减少地震数据和模型仿真数据之间的失配，其中在每次迭代中已更新模型的频谱被控制以匹配地下区域的输入的已知或估计频谱。
18.根据权利要求1所述的方法，其中所述成本函数是所述地震数据和所述模型仿真数据之间的交叉相关，且所述最优化在所述频谱控制应用到所述地震数据和所述模型仿真数据之后最大化所述成本函数。
19.根据权利要求18所述的方法，其中在所述地震数据上和所述模型仿真中使用源编码，且同步反演多个编码源。
20.根据权利要求19所述的方法，其中所述编码为所述反演的至少一次迭代变化。
21.根据权利要求1所述的方法，其中所述成本函数是所述地震数据和所述模型仿真数据之间的交叉相关的包络函数，且所述最优化最大化所述成本函数。
22.根据权利要求21所述的方法，其中在所述地震数据上和所述模型仿真中使用源编码，且同步反演多个编码源。
23.根据权利要求22所述的方法，其中所述编码为所述反演的至少一次迭代而改变。
【文档编号】G01V1/28GK103703391SQ201280016883
【公开日】2014年4月2日 申请日期:2012年1月30日 优先权日:2011年3月30日
【发明者】P·S·如斯, S·K·拉扎拉托斯, A·鲍姆斯泰因, I·希基舍韦, K·王 申请人:埃克森美孚上游研究公司