使用时变滤波器的全波场反演

摘要

本发明涉及一种在通过局部目标函数优化（64）执行地震数据（65）多尺度反演时，减小对起始模型的精度要求的改进方法。通过引入低通滤波器到目标函数（61），然后减小从一个尺度到另一个尺度被滤除的高频数据的量，实现不同尺度的反演。而且，滤波器被设计为是时变的，其中滤波器的低通截止频率随被滤除的地震数据的传播时间的增加而减小（62）。滤波器可如下设计：使用消除局部最小值的普拉特（Pratt）标准，仅针对源和接收器位置而非传播时间，执行周期和传播时间误差的平均（或其他统计度量）（63）。

说明书

相关申请的交叉参考

本申请要求2010年3月29日申请的美国临时专利申请61/318561 的权益，其标题为FULL WAVEFIELD INVERSION USING TIME  VARYING FILTERS，其全部内容通过引用包括在此。

技术领域

本发明一般涉及地震数据的数值转换从而推导传播介质的弹性参 数的领域。更具体地，本发明是用于在反演，如地震数据反演中执行 局部目标函数优化时降低对起始模型的精度要求的方法。

背景技术

反演（例如参考Tarantola，1984）试图发现最架解释观察的数据 的模型。最小化测量模拟数据和观察数据之间差的目标函数值的局部 反演方法通常是解决具有大量自由参数的模型的反演问题的唯一实用 方法。这些局部方法要求针对要反演的模型的初始猜测。这些局部方 法迭代更新模型，从而通过在基于目标函数梯度的方向上搜索当前模 型的扰动使其更接近真解。遗憾的是，目标函数通常具有许多最小值， 而不仅仅是一个对应于解模型的最小值。这些其他最小值被称为局部 最小值，而对应于所需解的最小值被称为全局最小值。如果反演的起 始模型太接近对应于一个这些局部最小值的模型，则局部反演方法在 局部最小值附近卡壳，并且永远也不能从该值迭代到全局最小值。因 此，无论如何努力迭代，都会产生错误解。

该局部最小值问题可通过对改变的目标函数执行第一迭代反演解 决，改变的目标函数具有较少的局部最小值但在所需解的位置附近有 全局最小值。对改变的目标函数迭代的结果应该产生更接近所需解的 模型。然后该更精确模型被用作针对原始目标函数上的迭代的初始模 型。因为这个新初始模型接近原始目标函数的全局最小值，所以对原 始目标函数的迭代现在应产生精确解。在改变的目标函数上迭代的这 个技术通常被称为多分辨率，或多网格，或多尺度反演，这将在下面 进一步讨论。

存在大量已知的反演方法。这些方法落入两个类别，迭代反演和 非迭代反演。下面是这两类中每一个通常表示的定义：

·非迭代反演——通过采用某些简单背景模型和基于输入数据更 新模型实现的反演。该方法不使用更新的模型作为另一个反演步骤的 输入。对于地震数据的情况，这些方法通常被称为成像、迁移、衍射 层析（diffraction tomography）、线性反演或波恩（Born）反演。

·迭代反演——反演涉及重复改进地下特性模型，以便找到能满 意地解释观察的数据的模型。如果反演收敛，则最终模型更好地解释 观察的数据，并将更接近地逼近实际地下特性。迭代反演通常产生比 非迭代反演更精确的模型，但计算成本更高。

波反演意味着任何基于波模拟器的反演，如声学或地震反演。在 波反演中最常用的迭代方法是目标函数优化。相对模型M，目标函数 优化涉及目标函数S(M)的值的迭代最小化，该目标函数是计算的和观 察的数据之间失配的度量（这有时也称为成本函数）。计算的数据是 利用被编程为使用物理学决定的当前模型表示的介质中源信号传播的 计算机模拟的。该模拟计算可通过几种数值方法实现，包括但不限于 有限差、有限元或光线跟踪。在Tarantola的文献【Tarantola，1984】 之后，最常用的目标函数是最小二乘目标函数：

S(M)=(u(M)-d)^TC^-1(u(M)-d)，（1）

其中T表示矢量转置算符，且：

M=是N个参数的矢量[m₁,m₂,...,m_N]^T的模型，

d=测量的数据矢量（相对源、接收器和时间采样的），

u(M)=针对模型M的模拟数据矢量（相对源、接收器和时间采样 的），

C=协方差矩阵。

目标函数优化方法是局部的或全局的【Fallat等人，1999】。全局 方法仅涉及计算模型群体{M₁,M₂,M₃,...}的目标函数S(M)并从该模型 群体中选择近似最小化S(M)的包括一个或更多模型的集合。如果需要 进一步的改进，则新选择的模型集合可用作生成可再次相对目标函数 S(M)测试的新群体模型的基础。全局方法比局部方法更可能对正确的 解收敛，但应用于具有许多模型参数的多尺度反演问题成本太高。已 知的全局反演方法包括蒙特卡洛（Montre Carlo）、模拟退火、基因和 进化算法。

局部目标函数优化涉及：

算法1：用局部目标函数优化更新模型的算法。

局部反演方法比全局方法更有效，并因此是用于大尺度反演问题的唯 一实用方法。通用的局部目标函数反演方法包括最陡下降、共轭梯度 和牛顿方法。

应该注意，的计算，算法1第二步要求针对N个模型参数 m_i中每个参数计算S(M)的导数。当N非常大（约超过1000），如果 必须为每个模型参数执行该计算，则该计算极其耗时。幸运的是，伴 随法可用来一次为所有模型参数有效地执行该计算【Tarantola，1984】。 用于最小二乘目标函数和网格状模型参数化（gridded model  parameterization）的伴随法由下面的算法总结：

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

局部目标函数优化通常比全局目标函数优化便宜得多，但要求更 精确的起始模型。要求该更精确的起始模型，是因为目标函数通常具 有许多最小值，且局部优化方法通常发现这些最小值中最接近的一个。 对应于真实模型的最小值是所谓的全局最小值，且所有其他最小值被 称为局部最小值。如果起始模型不是最接近全局最小值，则局部优化 技术可能产生不精确的反演模型，其对应于最接近的局部最小值。这 在图1中示出，其中目标是针对具有两个参数m1和m2的模型M的 反演。虚周线110显示目标函数的值作为参数m1和m2的函数。全局 最小值120是由实线黑圆标记的，且两个局部最小值130和140由灰 色填充圆示出。反演从初始模型M⁽⁰⁾(150)开始，然后由局部优化迭代 一次模型M⁽¹⁾，直到模型M⁽³⁾(160)。无论再尝试多少次局部优化的迭 代，反演的模型都仅更接近M⁽³⁾附近的局部最小值130，而非接近全局 最小值120。

已经提出了几种试图克服该局部最小值问题的方法。如上所述， 许多这类方法在反演的早期迭代中都涉及对改变的目标函数的迭代。 选择这个改变的目标函数从而具有较少的局部最小值，但具有原始目 标函数的全局最小值附近的全局最小值。通过该方法，早期迭代将产 生这样的模型，虽然其不精确，但显著更接近原始目标函数的全局最 小值。图2示出对应于图1的局部优化，但使用具有较少局部最小值 的改变的目标函数。改变的目标函数具有全局最小值210（实黑圆）， 其靠近原始目标函数（交叉阴影圆）的全局最小值220，但与其不在相 同位置。从初始模型M⁽³⁾(230)开始，该模型与图1中使用的初始模型 相同，使用改变的目标函数的两次迭代导致模型M⁽²⁾(240)。该模型M⁽²⁾可用作初始模型，以便进一步迭代，但现在使用初始目标函数。这在 图3中示出，其中来自图2的迭代2模型（以灰色示出），再编号为 310并用作起始模型。现在迭代对靠近全局最小值220而非靠近图1中 局部最小值的模型M⁽⁴⁾(320)收敛。因为起始模型比原始开始模型更精 确，所以反演现在迭代到正确解。

通常在改变原始目标函数时，改变的目标函数中局部最小值的数 目与原始函数的全局最小值和改变的目标函数的全局最小值之间的距 离逆相关。因此，从具有较少数目的局部最小值和最不精确局部最小 值的目标函数开始，然后演进到具有增加数目的局部最小值和增加精 度的全局最小值的目标函数，然后以对原始目标函数迭代结束，对改 变的目标函数序列进行迭代是有利的。对具有几个局部最小值的改变 的目标函数执行初始迭代的方法通常被称为多尺度或多网格方法，且 该技术的流程图在图4中示出。

过程从步骤410开始，选择原始目标函数的变化从而优化。该改 变的目标函数取决于要拟合420的数据，并在步骤430迭代，直到改 变的目标函数在步骤440被发现被充分最小化。（该值小于所选最大 值或满足另一个停止条件）。当发生该情形时，在步骤450确定是否 当前反演模型充分最小化原始目标函数。如果没有，则过程返回到步 骤410且要么选择新改变的目标函数或原始目标函数来优化。最终， 过程结束（460）。

在解决地震全波场反演（“FWI”）的局部最小值问题的文献中已 经提出了两个改变的目标函数：

·高阻滤波器：Bunks的文献（Bunks等人，1995年）描述了通 过对测量的数据和用于计算模拟的地震数据的源特征（source signature） 都应用时间不变高阻滤波器（有时被称为低通滤波器，意味着使低于 其截止频率的频率通过并摒弃截止频率以上频率的滤波器），改变最 小二乘目标函数。然后，随着反演进行迭代，这些滤波器的高阻，即 截止频率增加，其中没有滤波器应用于最终迭代（没有滤波器对应于 原始目标函数）。已知如何设计这样的滤波器；参看Press等人的 Numerical Recipes in FORTRAN，The Art of Scientifc Computing， Cambridge University Press（1992）。这也可从下面的资源获得，如 Seismic Un*x（参看：http:／／www cwp mines edu/cwpcodes/）。

·剥层法：Maharramov的文献（Maharramov等人，2007）教导了 反演的初始迭代应局部化在浅层，且随迭代进行而延伸到深度范围。 相应地，当仅浅深度被反演时，数据中仅较短时间被反演，因为浅模 型可仅预测较短时间部分的数据。

一般地，如果起始模型足够精确，从而预测任何传播模式的传播 时间到该模式的半个周期内，则FWI收敛到全局最小值。这也可称为 普拉特（Pratt）的标准。（“实际上，对于地震波形反演，这意味着许 多波形能量必须被预测在观察波形的半个波长内（通过初始模型）； 如果不是，则当预测的波形匹配观察波形内的错误循环时，则获得最 小失配模型。”——Pratt，1999）。

本发明是当执行局部目标函数优化时降低对起始模型精度要求的 改进方法。

发明内容

本发明方法适用于任何基于波模拟器的反演，如声学或地震反演。 在其一个方面中，本发明是用于改变目标函数的特定方法，对于起始 模型中的给定不精确性，该目标函数减少发现全局最小值所需的迭代 数目。减少迭代数目将相应地减少成本和计算时间。改变包括引入时 变滤波器到目标函数中。选择该滤波器以便测量的数据和计算的数据 之间传播时间差的某统计度量小于数据主周期的某个分数（通常四分 之一）的占优周期。这意味着滤波器是高阻滤波器。进一步选择该滤 波器以便该滤波器的高切断频率随传播时间的增加而减小，使得其为 时变滤波器。

参考图6的流程图，在一个实施例中，本发明是反演测量的地震 数据（65）从而推断地下区域的物理特性模型的方法，包括通过用目 标函数的局部最小值优化（64）执行测量的地震数据的迭代、多尺度 反演而连续更新模型，该目标函数计算模型模拟的地震数据和测量的 地震数据之间的失配（61），其中变化的低通滤波器被用于通过过滤 失配计算中测量的和模拟的地震数据（62），将目标函数连续地从一 个尺度变到另一个尺度，所述滤波器是时变的，其中滤波器的低通截 止频率随地震数据的传播时间改变，所述地震数据在反演的一个或更 多尺度上被滤波。滤波器可以用Pratt的标准设计，用于减小局部最小 值的数目，针对源和接收器位置而非传播时间修改从而包括数据的统 计度量（如平均值），使得滤波器是时变的（63）。多尺度反演的最 后阶段优选使用未变的目标函数，且因此其更有效地被提供更精确的 起始模型以帮助其收敛到全局最小值，导致优化的物理特性模型（66）。

如同任何数据反演一样，该过程在实际应用中是高度自动化的， 即，借助根据本公开内容编程的计算机执行。

附图说明

本发明及其优势可通过参考下面的详细描述和附图更好地被理 解，其中：

图1是收敛到局部最小值的反演的示意图；

图2是对应于图1的局部优化的示意图，但其使用具有较少局部 最小值的已改变的目标函数；

图3示出利用来自图1的原始目标函数的局部优化，其使用来自 图2的迭代2模型（灰色示出）作为起始模型；

图4是示出多尺度优化中的基本步骤的流程图；

图5是图解说明测量的数据和模拟的数据之间传播时间误差t_error， 和地震道的瞬时周期T的示意图；以及

图6是示出本发明一个实施例中基本步骤的流程图。

下面结合示例性实施例描述本发明。然而，一定程度上，下面的 详细描述是针对本发明特定实施例或特殊使用的，这仅是为了说明的 目的，而不能解读为对本发明范围的限制。相反，本发明涵盖权利要 求限定的本发明范围内的所有替换、修改和等价物。

具体实施方式

数学上，Pratt的标准可表述为：

$\frac{\underset{s,r,t}{\max }\left|t_{\mathrm{error}}(M,s,r,t)\right|}{\underset{s,r,t}{\max T(s,r,t)}}\le \frac{1}{2},---\left(2\right)$

其中t_error是测量的数据和模拟的数据之间的传播时间误差，且T是测 量的数据的瞬时周期，如图5中所示。在图5中，测量的数据和模拟 的数据之间的传播时间误差被指示为t_error，且地震道的瞬时周期由T 指示。传播时间误差是对齐测量的数据和模拟的数据所需的时移量。 瞬时周期是数据的相似相位（如，峰或槽）之间的时间。（传播时间 意味着从地震源发生地震波开始直到地震波在接收器被记录逝去的时 间。）

传播时间误差和瞬时周期都是源s、接收器r和传播时间t的函数。 此外，t_error是当前模型M的精度的函数。实际上，等式2可比必需的 约束更大，以确保FWI的收敛。特别地，可使用比最大值（如，平均 值）不严格的统计度量，或非1/2的目标可用在不等式的右侧。因此， 实际上，等式2可由下式取代：

$\frac{\underset{s,r,t}{\max }\left|t_{\mathrm{error}}(M,s,r,t)\right|}{\underset{s,r,t}{\mathrm{statT}(s,r,t)}}\le \frac{\alpha }{2}---\left(3\right)$

其中stat是某个统计量，如平均值、众数或均方差，且α是近似等于 一的常数。结果对统计的选择非常灵敏是没有预计到的。

通过利用高阻滤波器以增加T(s,r,t)，因此允许t_error的较大值，目标 函数的Bunk的变化遵从类似于Pratt的标准的推理逻辑。可以理解， 局部最小值的主要原因是跳周（cycle skipping），且较长周期导致这个 可能性较低。理论上，t_error可减小，而非限制数据到较低频率；然而， 实现这一点的唯一方式是具有更精确的起始模型，这非常困难并可能 是不可能的。而且，FWI的目标是产生更精确的模型，因此要求非常 精确的起始模型减小了FWI的值。另一方面，Maharramov的剥层法通 过仅反演浅模型来避免大的传播时间误差，该浅模型仅传播具有少的 传播时间的模式。因为传播时间误差通常随传播时间增加，所以限制 反演到较短传播时间将t_error保持在拇指规则内。

在本发明中，提出第等式3的“拇指规则”的替换，这导致确保收 敛的新策略。该拇指规则的替换规则如下：

$\frac{\underset{s,r}{\max }\left|t_{\mathrm{error}}(M,s,r,t)\right|}{\underset{s,r}{\mathrm{statT}(s,r,t)}}\le \frac{\alpha }{2}---\left(4\right)$

该拇指规则不同于等式3，因为统计分析不再对时间执行。在应用统计 计算后，等式4左侧的分子和分母不是源位置和接收器位置的函数。 等式4等价于：

$\hat{T}\left(t\right)\ge \frac{2}{\alpha }{\hat{t}}_{\mathrm{error}}(M,t)---\left(5\right)$

其中$\hat{T}\left(t\right)=\underset{s,r}{\mathrm{stat}}T(s,r,t),$且${\hat{t}}_{\mathrm{error}}(M,t)=\underset{s,t}{\mathrm{stat}}\left|t_{\mathrm{error}}(M,s,r,t)\right|.$

实际上，是传播时间t的增函数，因为传播时间误差倾向于随波 传播通过不精确的模型而累积。等式5表明针对多尺度反演的最优策略 会改变测量的数据和模拟的数据，以便地震数据的平均瞬时周期以类似 于的方式随传播时间增加。这可通过对测量的数据应用时变高阻滤波 器实现。这个时变滤波器的高截止频率应随增加的传播时间而减小。该 方法相对Bunk技术的优点是更多信息（即，在少的传播时间的较高频率 信息）会用于反演的早期迭代，因此更好地约束反演的模型，导致更快 的收敛。我们的建议相对Maharramov的剥层的优点是更多数据（即，传 播通过模型较深部分的地震模式）用于早期迭代，导致更快的收敛。

函数取决于传播时间误差是如何测量的和什么统计度量被应用 于这些传播时间误差测量两者。然而，可以预期，这个提出的多尺度反 演测量将对不够灵敏。实际上，不是努力从数据测量替换策略 是为假设简单的函数形式，如线性函数其中M₀是初始 模型。这个假定的函数形式然后被用于设计满足等式5的时变高阻滤波 器，且尝试使用该滤波器的反演。如果反演不收敛，则增加β的值，并 可尝试具有更保守估计的反演。继续这个测试直到发现为初始模型产 生收敛的β。

在发现对于这个初始模型，M₀收敛的β之后，迭代然后将产生比初 始模型更精确的当前反演的模型。这个增加的精度表明β应该随迭代的 进行减小。β的这个减小意味着通过更高频率的相应时变滤波器。反演利 用通过越来越多高频率的时变滤波器进行，直到数据中所有频率都用于 反演的最终迭代中。

为了使反演中的时变滤波器实用，能够以与时变滤波器一致的方式 计算伴随梯度（算法2）是必要的。完成这一点的数学上最直接方式是用 等式1中的逆协方差矩阵C^-1实现时变滤波器。为了完成这一点，选择逆 协方差矩阵C^-1为非对角线的（在时间维中），其中元素等于时变滤波器 系数的时间表示。因为滤波器会是时变的，所以滤波器系数随时间改变。 例子1示出示例性的C^-1的子矩阵，其对应于实现时变滤波器的特殊源和 接收器。该子矩阵的第一行为零，预期对角线上为一个1。

例子1

这意味着该特殊时变滤波器在时间零不执行滤波。在子矩阵中，非 对角线元素随行增加而增加，这意味着该时变滤波器的高截止频率随时 间增加而减小。注意，传播时间误差可被视作源或接收器的函数。在该 情形中，该方法可以以比简单时变滤波器更一般的方式被应用。例如：

1.传播时间误差通常不仅是传播时间的函数，而且通常是源和接收 器之间偏移（offset）的函数。如果是该情形，则为不同源接收器偏移使 用不同时变滤波器是有利的。随着迭代的进行，滤波器由于模型精度增 加而得到弛豫（relaxed）。

2.传播时间误差通常是源和/或接收器域中数据的时间倾斜（time  dip）的函数。这会发生是因为陡峭的时间倾斜对应于主要在水平方向上 传播的波，这些波对初始模型的精度更灵敏。在该情形中，时变倾斜滤 波器（如，除去具有高时间倾斜的地震事件的频率-波数滤波器）可取代 时变时间滤波器使用。随着迭代的进行，倾斜滤波器将随模型精度增加 而弛豫。

在任何情形中，滤波器（如空间变化和时变滤波器、时变倾斜滤波 器等）可在协方差矩阵C^-1中实现，如上所述的，且然后如上所述继续梯 度计算。

优选方法可能是如果可为当前模型M估计则时变滤波器 可被设计成与等式5一致。否则，使用β的初始猜值，估计是 线性函数βt是合理的。再次，时变滤波器应设计被成与β的估值和等式 5都一致。如果这个反演收敛，则过程完成。如果该反演不收敛，则增加 β，且尝试另一反演。增加β的这个过程继续，直到实现收敛。

实际上，表示时变滤波器的矩阵C^-1是非常稀疏的，且因此其应用于 算法2的步骤3中的数据残余是通过应用时变滤波器算子最佳实现的， 而非通过矩阵乘法。实际上，这个反演方法不可能对用于实现时变滤波 器的方法非常灵敏。因此，时变滤波器可最有效地被实现为分窗口时间 不变滤波器。换句话说，数据将被分隔成时间窗口，然后不同的时间不 变滤波器将被应用于不同窗口。

上面的描述针对本发明特定实施例，这是为了说明本发明。然而， 对本领域技术人员来说，显然可以对这里所述实施例做出许多修改和变 化。例如，本发明方法不限于地震数据，并可应用于其中多尺度反演被 用来避免局部最小值的任何数据。所有这类修改和变化都在权利要求限 定的本发明的范围内。

参考文献

Bunks，C．，F．M．Salcck，S．Zaleski，G．Chavcn，1995，“Multiscalc seismic waveform  inversion,”Geophysics，60，pp．1457-1473．

Maharramov，M．,U．Albcrtin，2007，“Localizcd image-difference wave-equation  tomography,”SEG Annual Mccting Expanded Abstracts，San Antonio,2007，pp.3009-3013． Fallat，M.R.，Dosso，S.F.，“Gcoacoustic inversion via local，global，and hybrid algorithms,” Journalofthe AcousticalSocietyofAmerica105，pp．219-3230(1999)．

Pratt,R.G.,“Seismic waveform inversion in thc frcquency domain，Part 1：Theory and  vcrification in a physical scale model,”Geophysics64，pp.888-901(1999).

Tarantola,A.,”Inversion of seismic reflection data in the acoustic approximation,″ Geophysics 49，pp．1259-1266(1984)．