[0002] 本申请要求2010年12月1日提交的
发明名称为“Simultaneous Source Inversion for Marine Streamer Data with Cross-Correlation Objective Function(利用互相关目标函数进行的海洋拖缆数据同时源反演)”的美国临时
专利申请61/418,694和2011年7月20日提交的相同发明名称的美国临时专利申请61/509,904的权益，二者的全部内容均被引入本文作为参考。
发明领域
[0003] 本发明总体涉及地球物理勘探领域，更具体地涉及地球物理
数据处理。具体地，本发明是得自多种地球物理源如
震源的数据反演方法，包括地球物理模拟——在进行一次模拟中计算来自多个同时活动的地球物理源的数据。
[0004] 发明背景
[0005] 即使应用现代的计算能
力，
地震全波场反演仍是计算昂贵的工作。但是，用此方法获得详细的地下表示的益处预期超过该阻碍。导致更快周转时间的
算法和工作流程的发展是使此技术适于场大规模数据的关键步骤。地震全
波形反演包括数据的正向和伴随模拟的数次
迭代(重复，iteration)。因此，降低正向和伴随计算运行成本的技术将致使用户在合理的时间量内解决较大规模的问题。
[0006] 地球物理反演[1，2]尝试找到最优解释观测数据并满足地质和地球物理限制的地下性质的模型。存在大量众所周知的地球物理反演方法。这些众所周知的方法属于两种类别即迭代反演和非迭代反演中的一种。以下是两种类别中每一种的普遍意义的定义：
[0007] 非迭代反演——通过假设一些简单的背景模型和基于输入数据更新模型实现的反演。此方法不将更新模型用作另一反演步骤的输入。对于地震数据的情况，这些方法常被称为成象、偏移、绕射层析成象或博恩反演。
[0008] 迭代反演——涉及重复改进地下性质模型从而建立理想地解释观测数据的模型的反演。如果反演收敛，则最终模型将更好地解释观测数据并将更紧密地接近实际的地下性质。迭代反演通常生成比非迭代反演更准确的模型，然而计算起来昂贵得多。
[0009] 迭代反演总体上相对于非迭代反演是优选的，因为其生成更准确的地下参数模型。不幸地，迭代反演的计算
费用如此昂贵，以致于将其应用于多种感兴趣的问题是不切实际的。这种高计算成本是源于所有反演技术均需要多个计算密集型模拟的事实。任何单独模拟的计算时间与所要反演的源数成比例，并且地球物理数据一般存在大量源，在此前述使用的术语源是指源设备的激活
位置。迭代反演的问题加剧，因为必须计算的模拟数与反演的迭代数成比例，并且所需迭代数一般属于成百上千的级别。
[0010] 用于地球物理学的最常用迭代反演方法是成本函数优化。成本函数优化包括成本函数S(M)的数值相对于模型M的迭代最小化或最大化，该成本函数S(M)是计算数据与观测数据之间错配度的度量(其有时也被称为目标函数)，其中计算数据是通过计算机使用当前地球物理性质模型和源
信号在由给定的地球物理性质模型表示的介质中的物理学支配传播(physics governing propagation)来模拟的。模拟计算可通过几种数值方法中任一种完成，该数值方法包括但不限于有限差分、有限元或射线追踪。模拟计算可在
频率域或时间域中进行。
[0011] 成本函数优化方法是局部或全局的[3]。全局方法简单地包括计算模型群体{M1，M2，M3,…}的成本函数S(M)和从大约最小化S(M)的群体选择具有一个或多个模型的组。如需进一步改进，则此新
选定的模型组可被用作生成新模型群体的
基础，该新模型群体可关于成本函数S(M)再次进行测试。对于全局方法，测试群体中的各模型可被认为是迭代，或在较高
水平上各组所测群体可被认为是迭代。公知的全局反演方法包括蒙特卡罗(Monte Carlo)、模拟
退火、遗传和
进化算法。
[0012] 不幸地，全局优化方法一般收敛极其缓慢，因此大多数地球物理反演是基于局部成本函数优化。算法1概括了局部成本函数优化。
[0013]
[0014] 算法1–进行局部成本函数优化的算法。.
[0015] 此程序通过利用新更新的模型作为另一梯度搜索的起始模型进行迭代。该过程持续，直到发现理想地解释观测数据的更新模型。常用的局部成本函数反演方法包括梯度搜索、共轭梯度和
牛顿法。
[0016] 在声学近似中地震数据的局部成本函数优化是普遍的地球物理反演作业，并通常是其他类型的地球物理反演的示例。当在声学近似中反演地震数据时，成本函数可写作：
[0017]
[0018] (方程式1)
[0019] 在此：
[0020] S＝成本函数，
[0021] M＝N个描述地下模型的参数(m1，m2，…mN)的向量，
[0022] g＝道集(gather)指数，
[0023] wg＝道集g的源函数，其是空间坐标和时间的函数，对于点源而言其是空间坐标的δ函数，
[0024] Ng＝道集数，
[0025] r＝道集中的接收器指数，
[0026] Nr＝道集中的接收器数，
[0027] t＝道(trace)中的时间样本指数，
[0028] Nt＝时间样本数，
[0029] W＝最小化标准函数(我们通常选择W(x)＝x2，其是最小二乘方(L2)标准)，[0030] ψcalc＝从模型M计算的地震压力数据，
[0031] ψobs＝测量地震压力数据。
[0032] 该道集可以是可在一轮地震正向建模程序中模拟的任何类型的道集。通常道集相应于地震炮点(shot)，尽管炮点可比点源更普遍(general)。关于点源，道集指数g相应于各个点源的位置。关于平面波源，g将相应于不同的平面波传播方向。这种广义源(generalized source)数据ψobs可以在现场得到或可由利用点源得到的数据合成。另一方面，计算数据ψcalc通常可通过利用正向建模时广义源函数直接计算。关于多种类型的正向建模——包括有限差分建模，广义源所需的计算时间粗略地等同于点源所需的计算时间。
[0033] 方程式(1)可简化为：
[0034]
[0035] 其中，现暗指基于接收器和时间样本的总和，并且
[0036]
[0037] 反演试图更新模型M，以使S(M)是最小的。其可通过局部成本函数优化实现，该局部成本函数优化如下更新给定模型M(k)：
[0038]
[0039] 其中，k是迭代数，α是模型更新的标量大小， 是关于模型参数取定的错配度函数的梯度。模型扰动或模型更新值通过目标函数梯度乘以步长α计算，其必须被反复计算。
[0040] 基于方程式(2)，可得到下列成本函数梯度方程式：
[0041]
[0042] 因此，为计算成本函数的梯度，必须分别计算各道集对成本函数的贡献梯度，然后加和那些贡献。因此，计算 所需的计算量(computational effort)是确定单个道集对梯度的贡献所需的计算量的Ng倍。对于地球物理学问题，Ng通常相应于地球物理源数，并且属于10,000至100,000的级别，大幅扩大了计算 的成本。
[0043] 注意， 的计算需要计算关于N个模型参数中的每一个mi的导数W(δ)。由于对于地球物理学问题而言N通常很大(通常大于一百万)，这种计算可能极其耗时——如果其必须对于每一个单独的模型参数进行。幸运地，可利用伴随(adjoint)法同时(at once)对所有模型参数有效地进行这种计算[1]。通过如下算法总结用于最小二乘方目标函数和
网格模型参数化的伴随法：
[0044]
[0045]
[0046] 算法2–利用伴随法计算网格模型的最小二乘方成本函数梯度的算法。
[0047] 虽然利用伴随法计算梯度相对于其他方法有效，但是其仍十分昂贵。具体地，伴随法需要两种模拟，一种在时间上正向，一种在时间上反向，而对于地球物理学问题这些模拟通常是极其计算密集型的。而且，如上所述，这种伴随法计算必须分别对每一个测量数据道集进行，增加了Ng倍的计算成本。
[0048] 所有反演种类的计算成本可通过自源组合反演数据而非单独地反演源而降低。这可被称为同时源反演。几种类型的源组合是已知的，包括：相干(coherently)加和近间距的源以生成有效源，该有效源生成一些所需形状(例如，平面波)的波前、加和宽间距的源、或完全或部分地
叠加反演前的数据。
[0049] 通过反演组合源实现的计算成本降低由于组合数据的反演通常产生不太准确的反演模型的事实而至少部分抵消。这种准确性的损失是缘于如下事实：当加和各个源时丢失信息，因此加和数据限制反演模型的强度不如未加和的数据。加和过程中的这种信息丢失通过在加和前编码各炮点记录被最小化。组合前进行编码显著保留同时源数据中的更多信息，因此更好地限制反演[4]。编码还允许组合近间距的源，因此允许给定计算区域的更多源组合。多种编码方案可用于这种技术，包括
时移编码和随机
相位编码。此背景部分的其余内容简要地评述多种已公开的地球物理同时源技术——编码和非编码的。
[0050] Van Manen[6]提出利用地震干涉测量方法
加速正向模拟。地震干涉测量通过将源布置在目标区域的边界各处而工作。对这些源单独建模，并且记录需要格林函数(Green’s function)的所有位置的波场。任两个记录位置之间的格林函数然后可通过互相关在两个记录位置处得到的道并在所有边界源基础上加和而计算。如果待反演的数据在目标区域内具有很多的源和接收器(相对于在边界上具有一个或另一个)，则这是计算所需格林函数十分有效的方法。但是，对于地震数据的情况，待反演数据的源和接收器均在目标区域中是很少的。因此，这种改进对于地震反演问题具有非常有限的应用性。
[0051] Berkhout[7]和Zhang[8]提出反演总体上可通过反演非编码同时源而改进，该非编码同时源在地下一些区域中被相干加和以生成一定程度上所需的波前。例如，点源数据可随时移而加和，该时移是源位置的线性函数，从而以相对于地表的一定特定
角度生成下行平面波。这种技术可用于所有种类的反演。此方法的问题是，源道集的相干加和必定降低数据的信息量。因此，例如，加和生成平面波去除了地震数据中所有关于传播时间相对于源-接收器炮检距(source-receiver offset)的信息。该信息对于更新缓慢改变的背景速度模型是至关重要的，因此Berkhout法没有很好地进行限制。为克服这个问题，可反演数据的多个不同相干加和(例如，多个具有不同传播方向的平面波)，但是然后由于反演成本与反演的不同加和数量成比例而失去效力。在本文中，这种相干加和的源被称为广义源。因此，广义源可以是点源或生成具有一定所需形状的波前的点源总和。
[0052] Van Riel[9]提出通过如下进行的反演：输入地震数据的非编码叠加或部分叠加(关于源-接收器炮检距)，然后限定关于此叠加数据的成本函数，其将被优化。因此，此公开提出利用非编码的同时源改进基于成本函数的反演。与Berkhout[6]的同时源反演方法一样，通过此方法提出的叠加减少了待反演数据中的信息量，因此该反演比对于原始数据得到较不充分的限制。
[0053] Mora[10]提出反演数据，该数据是宽间距的源的总和。因此，此公开提出利用非编码的同时源模拟提高反演效力。加和宽间距的源具有这样的优势：比Berkhout提出的相干加和保留更多信息。但是，宽间距的源的加和意味着必须用于反演的孔隙(aperture)(反演的模型区域)必须增加以容纳所有宽间距的源。由于计算时间与此孔隙的面积成比例，Mora法没有产生如加和的源彼此接近时所可以实现的一样多的效力增加。
[0054] Ober[11]提出通过利用同时编码源加速地震偏移——非迭代反演的特殊情况。在测试不同编码方法后，Ober发现，所得偏移图像具有显著降低的
信噪比，这是因为宽带编码函数必须仅大约
正交的事实。因此，当加和大于16个炮点时，反演
质量不理想。由于非迭代反演开始不是非常昂贵并且由于需要高信噪比反演，因此该技术在地球物理学工业中未被广泛实践。
[0055] Ikelle[12]提出通过同时模拟在不同时间间隔(在模拟中)激活的点源而快速正向模拟的方法。还讨论了将这些时移的同时源模拟数据解码返回单独的模拟的方法，该单独的模拟将已得自各个点源。然后这些解码数据可被用作任何常规反演程序的部分。Ikelle方法的问题是，提出的解码方法将生成单独的数据，其噪声水平与相邻源的数据之间的差异成比例。此噪声对于非横向恒定的地下模型将是显著的，例如包含倾斜
反射器(dipping reflectors)的模型。此外，此噪声将与同时源数成比例增长。基于这些困难，如果用于反演非横向恒定的地下，Ikelle的同时源方法可导致不可接受水平的噪声。
[0056] Krebs等在PCT专利申请公开号WO 2008/042081——本文将其被赋予的全部权限引入作为参考——中提出的源编码，是反演全波场数据非常成本有效的方法。(同时反演编码的道集的同样的方法将服务于接收器——通过源-接收器互易性或通过编码共源数据道集中的实际接收器位置)。对于固定接收器，单个有效源仅需进行正向和伴随计算；参见PCT专利申请公开号WO 2009/117174，本文将其被赋予的全部权限引入作为参考。如果事实是针对一般的2D采集几何记录数百炮点和在3D考察的情况下记录数千炮点，此方法带来的计算保存(saving)相当大。在实践中，对于最常见的场
数据采集几何，固定接收器假设不严格地有效。在海洋拖缆数据的情况下，每一个新炮点的源和接收器均发生移动。即使是在接收器位置固定的考察中，实践通常是：不是所有接收器都“监听”每一个炮点，并且正在监听的接收器可自炮点-至-炮点而不同。这也有悖于“固定接收器假设”。此外，由于逻辑问题，难以记录接近源的数据，并且这意味着一般缺少炮检距附近的数据。对于海洋和陆地考察也是这样的。这些因素均意味着，对于同时源道集，每一个接收器位置对于一些源炮点都将是漏测的数据。总之，在同时编码源反演中，对于给定的同时编码道集，需要每一个炮点的所有接收器位置处的数据，并且这可被称为同时编码源反演的固定接收器假设。在WO 08/042081中，当固定接收器假设不理想时，一些公开的实施方式可比其它的好。因此，对同时编码源(和/或接收器)反演的直接应用进行调节或调整将是有利的，这在固定接收器假设被折衷时将增强其性能。本发明提供了如此实施的方式。Haber等[15]也描述了关于利用随机优化方法移动同时编码源反演中的接收器的问题的方法，并将其应用于直流
电阻性问题。
[0057] 发明概述
[0058] 在一个实施方式中，本发明是利用同时源编码进行全波场反演以确定地下区域的物理性质模型的计算机执行方法，该全波场反演是来自不满足同时编码源反演固定接收器假设的考察的测量地球物理数据的全波场反演，该方法包括参考图2
流程图描述的下列步骤——其中所有加和、模拟、计算和更新在计算机上进行：
[0059] (a)获得测量地球物理数据(40)的两个或更多个编码道集组，该编码道集在固定接收器假设无效的情况下获得，其中各道集与单个广义源或可选地与单个接收器相关，并且其中各道集用选自非同等编码函数(30)组的不同编码函数编码(60)；
[0060] (b)通过加和相应于单个接收器(60)或可选地单个源的各道集中的全部或选定的数据记录，加和组中的编码道集，并针对每个不同接收器或可选地每个不同源重复进行，生成同时编码道集(80)；或可选地在场考察(90)中获得编码数据，其中编码函数是测量源特征(50)；
[0061] (c)假设地下区域的物理性质模型(10)，所述模型提供整个地下区域的位置处的至少一种物理性质的数值；
[0062] (d)利用假设的物理性质模型模拟数据的合成同时编码道集，其中模拟利用编码源特征(70)，并且其中在单个模拟操作中模拟整个同时编码道集；
[0063] (e)计算度量测量数据同时编码道集与模拟同时编码道集之间互相关的目标函数(100)；
[0064] (f)通过优化目标函数(100)更新(110)物理性质模型；
[0065] (g)利用自前次迭代的更新的物理性质模型(20)作为步骤(c)的假设物理性质模型，再重复步骤(a)-(f)至少一次，生成进一步更新的物理性质模型；和[0066] (h)在收敛(120)后，下载、显示或保存进一步更新的物理性质模型于计算机
存储器(130)。
[0067] 该方法进一步包括如步骤(a)所述获得至少一个另外的编码道集组，所述编码道集组具有所述测量地球物理数据的两个或更多个编码道集，和对于每个另外的组进行步骤(b)，然后累积来自步骤(e)的所述相应计算目标函数，其中步骤(f)所述的所述物理性质模型更新通过最大化所述累积计算目标函数而确定，其中测量数据的所述编码道集通过如下被编码：利用对于所述道集选择的所述编码函数在时间上褶积来自所述道集的所有道，其中用于所述合成同时编码道集模拟的所述编码源特征采用与用于编码测量数据的所述同时编码道集相同的编码函数，或是通过利用与用于编码步骤(a)的相应测量道集相同的编码函数在时间上褶积测量或估测的源特征而作成的函数，其中所述测量地球物理数据包括每个源激活的测量或估测的源特征，其中所述编码函数是选自线性、随机相位、线性调频脉冲、改良线性调频脉冲、随机时移和频率不相关的相位编码的类型，其中所述编码函数对于一些源是一种类型，并且对于其他源是另一种类型，其中所述编码函数被优化以提高所述目标函数的质量，其中步骤(d)所述的模拟利用有限差分、有限元或有限体积模拟代码进行，其中所述物理性质模型是
地震波速度、地震弹性参数、地震
各向异性参数或地震滞弹性参数的模型。
[0068] 在一个实施方式中，本发明是从地下区域生产
烃的方法，包括：
[0069] (a)进行所述地下区域的地震考察，其中不满足同时编码源反演的固定接收器假设；
[0070] (b)获得所述地下区域的速度模型，所述速度模型通过包括如下的方法确定：
[0071] 通过迭代反演，反演考察测量数据，所述迭代反演包括，利用假设或更新的速度模型同时模拟代表多个考察源或可选地多个接收器的考察数据，其中所述模拟中的源或接收器特征被编码，生成地球物理数据的模拟同时编码道集，所述反演进一步包括，通过优化度量所述模拟同时编码道集与所述测量数据相应的同时编码道集之间的互相关的目标函数，确定用于下一次迭代的速度模型更新；
[0072] (c)将井眼钻探到至少部分由利用(b)的更新速度模型完成的对所述地下区域结构的阐释而确定的所述地下区域中的层中；和
[0073] (d)从所述井眼生产烃。
[0074] 其中所述编码通过选定的编码函数将各道集和各源特征相乘而进行，并且其中在步骤(g)的每次迭代中，一个或多个编码函数被设定为等于零，以使不归零的所述数据记录相应于平均比编码前的所述测量地球物理数据相距更远的源激活位置，并且随着所述迭代进行所有数据记录均被应用，即，不归零。
[0076] 由于专利法限制，一个或多个附图是彩色原图的黑白再现。彩色原图已提交于对应的美国申请中。具有彩色附图的此专利或专利申请公开副本可通过
请求和支付必要的费用从美国专利商标局获得。
[0077] 本发明及其优势将通过参考下文详细描述和附图而得到更好的理解，在附图中：
[0078] 图1A-1C是测试实例的数据，其中图1A显示应用中的“真实”速度模型，即，用于模拟合成考察数据的模型；图1B显示通过利用常规L2目标函数反演模拟考察数据获得的速度模型；和图1C显示通过利用本发明方法的互相关目标函数反演而获得的速度模型；
[0079] 图2是显示本发明方法的一个示例性实施方式的基本步骤的流程图。
[0080] 图3是显示间距20m的四个炮点的总和的模拟地震数据的截面；
[0081] 图4是显示间距100m的四个炮点的总和的模拟地震数据的截面；和[0082] 图5A-5C显示模型研究的结果，显示本发明方法的实施方式利用互相关目标函数减少在迭代反演中影响收敛的局部极小化问题的优势。
[0083] 本发明将结合实例实施方式进行描述。但是，在下文详述针对于本发明的具体实施方式或具体应用方面，这仅意为示例性，并不被解释为限制本发明的范围。相反，意为
覆盖可包括在所附
权利要求限定的本发明范围内的所有替代选择、修正和等同形式。
[0084] 实例实施方式详述
[0085] 本专利文本中提出的方法利用互相关目标函数反演编码同时源数据。经编码形成同时源数据的原始炮点道集具有漏测的近-和远-炮检距，这在海洋拖缆考察中是常见的。对于这种互相关目标函数如何用于非固定接收器数据的全波场反演(“FWI”)的直观理解可通过考虑成像问题而获得。在FWI中，反演一种或多种模式的地震波场，获得地球的物理性质模型。首先考虑同时源逆时偏移(“RTM”)的成像问题。
[0086] 假设源于两个源的两个地下测量波场d1,d2，并且其通过如下编码，[0087]
[0088] 其中P是将数据投影到固定接收器几何所需的所有接收器上的投影算子。因此，如果所有炮点的算子P相同，则是固定接收器几何。编码源由如下给出，
[0089] S＝α1S1+α2S2
[0090] 将正向问题表示为Ld＝S，其中L是偏微分算子，d是模拟响应，和S是编码源函数。如果将正向算子表示为L-1使得模拟数据d＝L-1S并将反向传播算子表示为 则通过互相关获得的图像由如下给出，
[0091]
[0092] 其中P是使全部区域的计算数据映射至固定接收器组的投影算子。对于固定接收器情况，所有源的P相同，因为每个接收器监听所有源。
[0093] 对于非固定接收器几何，假设Pk是第k个源到接收器子组上投影算子。则编码数据由如下给出，
[0094]
[0095] 则图像由如下给出
[0096]
[0097] 如果采用几种编码实现方式，则第二项(second term)由于串音而消除并且其对图像的贡献减少，因而增强信噪比。成像与生成FWI梯度之间存在紧密的联系。如果考虑FWI的目标函数(通常被称为成本函数)为L2，即，最小二乘方，测量数据与模拟数据之间的范数(norm)，然后L2目标函数的梯度，简单地是反向传播的数据剩余与模拟波场的互相关。对于固定接收器几何，这是良好的度量，因为特定接收器监听全部源。对于非固定接收器，这存在问题，因为编码测量数据相对于编码模拟数据是不足的。这种不足最终映射成模型假象。
[0098] 但是，从接收器非固定的成像条件获得的认识可被认为是使用互相关目标函数进行FWI可有助于减轻此问题的线索。因此，替代L2，考虑互相关目标函数，其由如下给出：
[0099]
[0100] 其中 是非固定接收器的编码测量数据，暗示将接收器未监听源的道假设是零道。 是来自同时源函数的模拟编码数据；这是固定接收器数据，因为所有源同时被爆破(fired)，并且各接收器测量来自此同时源的信号。假设对于各炮点，假想地，我们能够分开考察中存在的接收器和漏测的接收器，从而可写成
[0101]
[0102] 贡献于 的道即是测量数据中不可用的道，即， 因此上述方程式中第二项的贡献相对于第一项很小。可将目标函数近似为，
[0103]
[0104] 而且，如果不同的编码用于一些迭代，优选在每次迭代中，如Krebs等在PCT专利申请公开号WO 2008/042081中教导，第二项最终减少。除每次迭代改变编码外，进一步的选择是在每次迭代应用编码组，该编码组被设计以获得目标函数的平均梯度，即，得自各编码的梯度在加和时给出平均梯度。因此，在每次迭代时，这种平均处理使串音减少——非常类似于同时源RTM成像的思想，在同时源RTM成像中，利用不同编码获得的图像叠加以消除串音噪声。梯度平均也在随机优化方法中被主张(Haber等，[15])。
[0105] 本发明方法的优选实施方式的总体目标是最大化——而非最小化——互相关目标函数。但是，可改写目标函数，例如，用其负数(negative)将其替代，在这样的方式下有效解决最小化问题。因此，本发明方法中此步骤的总体表述是优化目标函数。根据
优化算法的选择，即，其解决最大化还是最小化算法，可选择目标函数，以使其保留对任何任意数据规模不敏感的互相关性的主要特征。互相关目标函数的优势是其使FWI对小波规模不敏感。有效地，波形的传播时间要素适于此目标函数。如果目标函数以相位是传播时间的函数的频率域撰写，则此可容易看到。.
[0106] 合成实例
[0107] 本发明方法的可行性通过利用图1A所示速度模型的合成实例得以证明。在该实例中，建立具有80m炮点间距的383个炮点，并且每20m布置接收器。为模拟非固定接收器几何，对于各炮点去除近和远炮检距，并采用正200m至3km的炮检距。如预期，近和远炮检距的漏测接收器数据有悖于固定接收器假设。那些炮点经编码形成测量编码数据。在非线性反演的每次迭代时，应用考察中全部连续炮点，并用不同的编码函将其编码8次，以形成8个同时源，从而计算平均梯度。取决于平均梯度中存在多少串音，可增加编码数量。此被执行用于L2目标函数和互相关目标函数。
[0108] 利用带有
时间窗的多
分辨率法(multi-resolution method)实施反演。在多分辨率法中，数据通过频率带窗化，从而稳定化反演问题。一般，较低频率带的数据先反演，并且得自低频反演的模型用作下一个较高频率带的起始模型。该技术用于避免收敛于错误的解决方案——通常被称为局部极小化解决方案。图1B显示利用L2目标函数获得的模型。在模型的浅部和深部，假象均明确可见。图1C的模型是应用互相关目标函数的结果。模型清楚显示相对于图1B模型的显著改进，强化了互相关度量明显更适于处理非固定接收器几何同时源数据的观察结果。
[0109] 应用互相关目标函数的反演在所增加的源记录间距过近时可遭遇局部极小化。进一步拉开炮点记录的间距解决了此问题，但使加入反演的信息量减少，因而降低最终反演模型的质量。解决此矛盾的方式是将所有可用的炮点记录分成几组，其中进一步拉开属于同组的炮点记录的间距。然后可在非线性反演的每次迭代中选择不同炮点组，直到所有组均已被使用。此过程可重复，直到收敛。这可进行的一种方式是利用使源特征与测量数据道集相乘的编码函数。所有不在某次迭代选定组中的数据记录的编码函数可被设定为等于零。事实上，可利用编码函数随机选择各组中的炮点，按照Boonyasiriwat and Schuster，2010表现的思想，这将导致炮点间距增加，其中，它们将这用于固定接收器传播和不同的目的。这可利用随机编码顺序[-1,0,+1]进行。这种增加炮点间距的技术可用于本发明方法，以减少非固定接收器的串音噪声。可设计其他实现较小炮点记录组的较大炮点间距的方法。
[0110] 下文是应用上述策略的实例。
[0111] 此方法的可行性通过合成实例得以证明。建立了383个具有20m炮点间距的炮点，并且每10m布置接收器。为模拟非固定接收器几何，各炮点采用正200m至3km的炮检距。那些炮点经编码形成测量编码数据。在非线性反演的每次迭代时，利用8个同时源计算平均梯度。图3显示4个近间距(在此为20m)炮点的总和，而图4显示进一步远离(在此为100m)的4个炮点的总和。基本思想是炮点间距增加使目标函数的串音要素减少；但是，过大间距导致数据空隙，其降低更新模型的能力。如上所述，这种情况可通过选择若干组形成编码炮点而减轻。各组具有明显隔离的炮点，但采用这些组中的几组填补了空间空隙，否则如果仅选择一组会保留该空间空隙。对于在此呈现的实例，每次迭代选择4组，并将2种实现方式用于各组，实现方式是用不同编码函数组进行同时源模拟。因此，每次迭代总计有8个编码同时炮点。为进行比较，每次迭代以8种实现方式用近间距编码炮点进行反演。图5A是此实例所用的真实模型。图5B显示在源近间距时同时源反演的结果。由于模型中存在串音，反演超过50次迭代也无法建立适当的更新。50次迭代的模型被作图于图5B。为进行比较，炮点明显隔离时50次迭代的模型被作图于图5C。反演过程中不存在收敛问题。图5C与5A之间比图5B与5A之间具有明显更好的一致性。
[0112] 上述申请涉及本发明的具体实施方式，其目的是将本发明示例。但对于本领域技术人员明显的是，对本文所述实施方式的多种改动和变型是可以的。所有这种改动和变型均意为处于如所附权利要求限定的本发明范围内。本领域技术人员容易理解，在本发明的优选实施方式中，本发明方法的至少一些步骤在计算机中进行，即，本发明是计算机执行的。在这种情况下，所得更新物理性质模型可被下载、显示或保存于计算机存储器。
[0113] 参考文献
[0114] 1.Tarantola,A.,"Inversion of seismic reflection data in the acoustic approximation,"Geophysics 49,1259-1266(1984).
[0115] 2.Sirgue,L.,and Pratt G."Efficient waveform inversion and imaging:A strategy for selecting temporal frequencies,"Geophysics 69,231-248(2004).[0116] 3.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).
[0117] 4.Hinkley,D.and Krebs,J.,“Gradient computation for simultaneous source inversion,”PCT Patent Application Publication No.WO 2009/117174.[0118] 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.[0119] 6.Van Manen,D.J.,Robertsson,J.O.A.,Curtis,A.,"Making wave by time reversal,"SEG International Exposition and 75th Annual Meeting Expanded Abstracts,1763-1766(2005).
[0120] 7.Berkhout,A.J.,"Areal shot record technology,"Journal of Seismic Exploration1,251-264(1992).
[0121] 8.Zhang,Y.,Sun,J.,Notfors,C.,Gray,S.H.,Cherris,L.,Young,J.,"Delayed-shot 3D depth migration,"Geophysics 70,E21-E28(2005).
[0122] 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).
[0123] 10.Mora,P.,"Nonlinear two-dimensional elastic inversion of multi-offset seismic data,"Geophysics 52,1211-1228(1987).
[0124] 11.Ober,C.C.,Romero,L.A.,Ghiglia,D.C.,"Method of Migrating Seismic Records,"U.S.Patent No.6,021,094(2000).
[0125] 12.Ikelle,L.T.,"Multi-shooting approach to seismic modeling and acquisition,"U.S.Patent No.6,327,537(2001).
[0126] 13.Romero,L.A.,Ghiglia,D.C.,Ober,C.C.,Morton,S.A.,"Phase encoding of shot records in prestack migration,"Geophysics 65,426-436(2000).[0127] 14.Jing X.,Finn,C.J.,Dickens,T.A.,Willen,D.E.,"Encoding multiple shot gathers in prestack migration,"SEG International Exposition and 70th Annual Meeting Expanded Abstracts,786-789(2000).
[0128] 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).
[0129] 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,”Geophysics 74-6,WCC177-WCC188(2009).[0130] 17.Boonyasiriwat,C.,and Schuster,G.“, 3D Multisource full-waveform inversion using dynamic random phase encoding,”SEG Expanded Abstracts 29,1044-1049(2010).
