CN102884447B - Q层析成像方法 - Google Patents

Abstract

本发明涉及一种通过执行基于射线（60）的形心频率迁移（50）Q层析成像，从地震数据（10）重构地下Q模型（110）的方法。地震源波形的振幅谱由频率（40）的加权频率指数函数近似，具有两个参数以便调整以拟合频率迁移数据，从而提供对各种不对称源振幅谱的拟合。箱约束可用于优化程序，以及用于速度层析成像的多指标有效集方法是用于实施箱约束（100）的优选技术。

Description

Q层析成像方法

相关申请的交叉引用

本申请要求于2010年5月5日提交、题为“Q层析成像方法（QTOMOGRAPHY METHOD）”的美国专利临时申请61/331694的优先权，该申请的全文通过引用合并于此。

技术领域

本发明一般涉及地球物理探测领域，更具体涉及地震数据处理。具体地，本发明涉及Q层析成像的技术领域。

背景技术

针对岩石属性的特性描述和适当的振幅随偏移距变化（AVO）分析，需要考虑地震衰减效应。在迁移过程中，需要地震衰减信息来补偿吸收效应，以提高迁移图像质量。因此，地震衰减的估算对于储油层的检测和监视是必要的。

地震衰减可以由品质因数Q量化描述。简单假设地震衰减是频率相关的，但品质因数Q是频率不相关的。这个假设在勘探地球物理应用的频率范围内是有效的。Q层析成像是Q估算的方法。这个方法尝试从地震数据重构地下2D或3D的Q模型。通常，Q层析成像算法被分为两个主要类别。一个类别是基于射线的层析成像（Quan和Harris，1997；Plessix，2006；Rossi等人，2007）。另一个类别是基于波方程的层析成像（Liao和McMechan，1996；Hicks和Pratt，2001；Pratt等人，2003；Watanabe等人，2004；Gao等人，2005）。基于波方程的层析成像在物理上更为精确但是计算昂贵并且对3D情况不适用。本发明属于基于射线的Q层析成像类别。

Q层析成像的一个主要问题是如何以最小的近似值和最大的灵活性建立Q模型与地震数据之间的联系。广泛使用的方法是基于Q与地震振幅衰减之间的关系。另一个方法使用地震形心频率降档来估算品质因数Q。后面的这个方法被认为更加稳健，这是因为这个方法与几何扩散效应和反射/传输损失无关。不过，常规形心频率迁移方法仅能使用高斯函数、矩形波函数或三角形函数来拟合源振幅谱，这会引入明显的误差，这是因为在多数情况下，源谱不能由这些函数近似。本发明包括加权频率指数函数，其被设计为拟合各种不对称源振幅谱，从而通过大大减少源振幅谱拟合误差，改善Q层析成像的精度。

在多数现有Q层析成像算法中，优化部分是基于无约束的优化方法或基于简单的非负约束的优化方法。因此，这些Q层析成像算法花费大量计算时间或产生许多人工产物和不合实际的Q模型（例如，负的Q值或极端低的Q值），尤其是当地震数据被噪音污染时。包括带有箱约束的高效优化算法的本发明能够改善重构Q模型的质量和可靠性。下面进行现有技术的更详细讨论。

地震衰减层析成像（Q层析成像）已被研究许多年并取得很大进展。基于射线Q层析成像算法的两个主要组成部分是1）用于构建Q层析成像的数学模型的简单但是精确的地震数据与Q模型之间的关系；2）解这个数学问题的可靠且稳健的优化算法。已经开发或提议用于构建这两个组成部分的许多技术。这些技术讨论如下。

建立地震数据与Q模型之间的联系

估算Q的最简单直接的方法是谱比方法（Spencer等人，1982；Tonn，1991），其中两个地震波形之间谱比的对数被计算为频率的函数，并且这个函数由线性频率函数近似，其斜率被视为累积的地震衰减并最终与沿波传播路径的Q值相关。理想地，这个方法假设几何扩散效应和反射/传输损失效应与频率无关，消除了几何扩散效应和反射/传输损失效应。在实际应用中，由于子波重叠、线性拟合的不确定性以及许多其他因素，这个方法是相对不可靠的。

Rickett（2006）提出利用时间-频率分析技术辅助的谱比层析成像扩展方法。这个方法被宣称对绝对比例不敏感并应用于利用垂直地震剖面（VSP）的Q剖面估算的应用中。在这个方法中，描述与频率无关的振幅变化的对数-振幅标量包含在未知数中，这大大增加了未知数的数量并降低算法的效率。而且，在利用表面地震波反射数据的2D/3DQ层析成像中，对数-振幅标量不仅是位置的函数，而且是射线的函数，这使过程非常复杂。

基于地震子波上升时间的变化与沿传播路径的1/Q剖面是线性相关的，Wu和Lees（1996）报告了利用地震学中的上升时间的地震衰减层析成像方法。遗憾的是，这个方法在勘探地球物理中是不实用的，这是因为子波不可避免地被噪音、散射效应、重叠等污染。

应当指出，地震子波振幅谱的形状几乎完全受品质因数Q影响，并且开发了用于Q估算的峰值频率变化方法（Zhang，2008）。这个方法是有吸引力的，但是在实际中，存在精确测量峰值频率变化的困难。而且，由于仅使用在单个频率的信息，所以Q估算的不确定性大。

更稳健的方法由Quan和Harris引入（1997），其中，地震波形的整个频段的信息被用于计算形心频率降档，接着通过简单封闭形式公式，使形心频率迁移与沿射线路径的Q剖面相关。这个方法本质上是完全不受几何扩散和反射/传输损失影响的。这个方法的局限在于源振幅谱必须是高斯、矩形波或三角形函数。众所周知的是地震振幅谱从来都不是矩形波或三角形函数。而且，其通常是不对称的并且与高斯函数有很大不同。若这个不对称振幅谱由高斯函数近似，则会在Q模型的重构中引入明显误差。因此，若在不丢失记录的地震数据的形心频率与沿射线路径的Q剖面之间关系的简单性质的情况下，存在能够用于更精确拟合各种地震频谱的函数，则这个稳健方法在实际Q层析成像应用中会更加精确。

用于Q层析成像的约束优化算法

当建立地震数据与Q模型之间的关系时，基于射线的Q层析成像问题可以由线性优化问题描述。在大多数现有Q层析成像算法中，这个线性优化问题通过利用Krylov子空间方法迭代求解，例如共轭梯度方法和LSQR方法，而不应用任何约束（Quan和Harris，1997；Plessix2006，Rossi等人，2007）。只要地震数据具有高信噪比（SNR），那么这些算法工作得很好。不过，地震数据从来都不是干净的；在处理真实的场数据时，这些非约束的优化算法总是产生某些负的Q值或极端小的正Q值，这些是物理上不现实的。而且，在某些情况下，Q值范围的先验信息是已知的。在这些情况下，这个信息需要通过箱约束包含在层析成像算法中，以便提供更加可靠的Q层析成像结果。

Rickett（2006）开发出具有约束的Q估算算法。但是他的算法采用非负约束，而不是箱约束，这意味着消除了负的Q值，但仍然存在极低的正Q值。Rickett（2006）报告了应用非负约束的两种方法。第一种方法是非线性变换方法。在这个非线性变换方法中，变量Q由e^y代替，并且求解y而不是1/Q。通过这样做，Q被迫为正，但是整个***会是非常非线性的。为了实现这个目标，由此产生的优化***利用高斯-牛顿途径求解，而这是很昂贵的。这个方法的另一个缺点是，在优化过程中当Q值很大时，基于梯度的优化算法将是停滞的或收敛非常缓慢，即Q值会保持在那里并且不再变化。在最恶劣情况下，若在开始模型中的y值是无限的，那么成本函数的梯度是0并且优化算法不执行搜索。应用Rickett（2006）报告的非负约束的第二个方法是通过平滑技术，强制衰减的单调递增特性。这个方法利用VSP数据对Q估算是有效起作用的，但在利用表面反射地震数据的2D Q层析成像中会失效。

本发明人知道，不存在带有箱约束的现有Q层析成像算法来实施在上边界和下边界规定范围内的估算Q值。不过，在某些其他地球物理应用中，例如速度层析成像中采用带有箱约束的优化算法。例如，Delbos等人（2006）开发出具有箱约束的地震反射层析成像算法。在他们的算法中，利用高斯—牛顿增强朗格拉日方法求解约束优化问题，而关联的朗格拉日问题——另一个约束优化问题，通过梯度投影方法、有效集方法以及共轭梯度方法的组合来求解。他们使用的有效集方法是常规的，并且由于该算法更新有效集，一次一个约束，因此是低效的（Bierlaire等人，1991；1984；Nocedal和Wright，1999）。当箱约束的数量巨大时，该算法的收敛速率可以非常慢。

在本发明中，在数学领域中的最近进展（其被称为多指标有效集方法（Morigi等人，2007））被用于执行带有箱约束的Q层析成像，相比于非约束的Q层析成像算法和利用常规有效集方法的约束算法，其明显改善在Q重构质量和算法效率方面的Q层析成像算法的性能。

发明内容

在一个实施例中，本发明是基于射线的形心频率迁移Q层析成像方法，其从在利用地震源的勘测中由接收器测量的地震数据重构1/Q的地下深度模型，其包括选择数学函数来近似地震源的振幅谱，以便计算由于地球衰减引起的谱的形心频率迁移，并且将所述形心频率迁移与由品质因数Q的倒数表示的衰减相关，以及通过迭代的线性优化求解Q或1/Q，其中所述优化具有箱约束，以将估算的Q值保持在由上边界和下边界规定的与位置相关的范围内。该约束优化可以由多指标有效集方法求解，该方法允许一次多网格指标的有效集更新，其中网格指标表示地下位置。选择的数学函数可以是频率的加权频率指数函数。无论在优化中是否使用箱约束，使用这个函数来近似地震源的振幅谱将是有利的。

除了其他优点之外，由本发明方法产生的重构的地下Q模型可以有利地用于地震成像中，以补偿振幅模糊、频率损失以及由粘声波覆盖层，例如天然气储集层导致的相位失真效应。通过在地震成像过程中包含由本发明提供的更精确Q模型，能够明显改善地质结构图像的质量。此外，由于Q对某些岩石和流体属性，例如流体饱和度与孔隙度很敏感，所以重构的Q模型在储集层特征描述的应用中是有益的。在这样的应用中，本发明成为碳氢化合物勘探、开发或生产的方法。

Q层析成像领域的技术人员应该意识到本发明方法的至少某些步骤优选在计算机上执行，根据本文所述教导编程，即本发明在大多数或全部实际应用中是计算机实施的。

附图说明

通过参考下列详细描述及其随附绘图，可以更好理解本发明及其优势，其中：

图1示出本发明至少某些实施例中设计用于拟合各种不对称源振幅谱的加权频率指数函数的带宽变化，针对特性频率参数的各种值和称为对称指标的第二参数的选择值在n=2时的情况；

图2示出加权频率指数函数在n=1，2，3，4，5而形心频率（n+1）f₀固定在30Hz时的曲线图；

图3示出在Q层析成像中实施本发明的基本步骤的流程图；

图4到8A-B属于合成示例应用，其中

图4示出源振幅谱；

图5A-B示出速度模型和射线路径（5A）以及真实的1/Q模型（5B）；

图6A示出利用具有高斯函数拟合和没有箱约束的常规形心频率迁移Q层析成像方法的重构的1/Q模型；以及图6B示出图6A的重构1/Q模型与图5B的真实1/Q模型之间的差异；

图7A示出利用带有本发明的加权频率指数函数拟合、但是没有箱约束的形心频率迁移Q层析成像方法的重构的1/Q模型；图7B示出图7A的重构1/Q模型与图5B的真实1/Q模型之间的差异；以及

图8A示出利用带有本发明的加权频率指数函数拟合和根据本发明的箱约束的形心频率迁移Q层析成像方法的重构的1/Q模型；图8B示出图8A的重构1/Q模型与图5B的真实1/Q模型之间的差异。

图5A-B到8A-B是由于专利法在绘图上的约束，彩色显示的黑白复制品。可以通过从美国专利和商标局请求对应美国专利申请的副本并支付必要的费用获得彩色图的副本。

将结合示例性实施例来描述本发明。不过在某种程度上，下列详细描述具体到本发明的特定实施例或特定用途，这仅仅是为了说明目的，而不应解释为对本发明范围的限制。相反，正如所附权利要求所定义的，意欲覆盖包含在本发明范围内的所有替换、修改和等同物。

具体实施方式

本发明包括从地震数据重构2D或3D地震品质因数（Q）模型的方法，称为Q层析成像的技术领域。

本发明在至少某些实施例中的主要特征如下。专门设计的函数，加权频率指数函数分析和拟合源子波振幅谱。通过调整加权频率指数函数的两个参数实施源子波振幅谱拟合。已接收地震波形的形心频率相对源子波形心频率的迁移被计算并输入到带有箱约束的优化算法中，以便重构Q模型，其中Q值的范围由先验信息预先确定。不像最初的形心频率迁移方法（Quan和Harris，1997），利用专门设计的振幅谱拟合函数，这个Q层析成像算法能够处理非高斯且非对称源振幅谱，以减少由源振幅谱拟合中的不匹配引入的误差。利用这个专门设计的源谱拟合函数，Q层析成像问题转换为带有箱约束的约束优化问题。这个约束优化通过采用多指标有效集方法（Morigi等人，2007）求解，其进一步改善这个Q层析成像算法的精度和稳健性，而不牺牲高效率特征。术语有效集指的是指示的不能够在迭代周期结束时更新的正被优化未知量组的子集，这是因为他们碰到约束，无论是上限还是下限。

下一步解释本发明的某些基本理论。

如果源子波的振幅谱是S(f)，那么接收的地震波形的振幅谱R(f)可以表述为（Quan和Harris，1997）

R(f)＝GH(f)S(f)            (1)

其中G是包括几何扩散效应、反射/传输系数等的、与频率无关的因数。H(f)是描述地震衰减效应的脉冲响应函数，其被公式化为

$H\left(f\right)=\exp (-f{\∫}_{\mathrm{ray}}\frac{\mathrm{\π}}{\mathrm{Qv}}\mathrm{dl})---\left(2\right)$

其中Q是品质因数以及v是地震（即，声学）波速度。

在用于Q估算的形心频率迁移方法中，关键部分是使用分析函数来拟合源振幅谱S(f)，接着导出已接收地震的形心频率与沿波传播路径的Q剖面之间的显性关系。在最初的形心频率迁移方法中，这个显性关系仅通过假设源振幅谱是高斯函数或矩形波函数（除了单个间隔处是等于常数以外，在整个实线上是零的函数）或三角形函数导出。这是传统形心频率迁移方法的主要不足，其会导致Q估算中的大量误差（Rickett，2006；Zhang，2008）。

本发明至少某些实施例的一个主要特征是加权频率指数函数被设计用于更精确拟合各种不对称源振幅谱，而不丢失形心频率迁移与沿射线路径的Q剖面之间的简单封闭形式关系。而且，利用这个专门设计的函数实施源幅度谱拟合是非常方便的，这是因为函数的形状和带宽由两个参数分别确定。加权频率指数函数的公式是

$F\left(f\right)=A{f}^n\exp (-\frac{f}{f_0})---\left(3\right)$

其中A是振幅比例的常数，f₀称为特征函数，以及n称为对称指标。特征频率f₀是用于带宽控制的参数。如果对称指标n固定，则这个函数的带宽随着特征频率增加而扩展。实际上，F(f)的形心频率具有非常简单的形式，给出如下：

$f_F=\frac{{\∫}_0^{+\∞}\mathrm{fF}\left(f\right)\mathrm{df}}{{\∫}_0^{+\∞}F\left(f\right)\mathrm{df}}=(n+1)f_0---\left(4\right)$

图1示出当的n=2时F(f)的带宽变化，具有范围从10Hz到50Hz的增量为10Hz的五个不同特征频率。在图1中，当f₀从50Hz变化到10Hz时，F(f)的形状保持相对不变，而带宽如预计的那样萎缩。另一方面，对称指标n被用于控制F(f)的对称属性。对称指标n越大，函数F(f)的形状越对称，如图2所示，其中形心频率固定为30Hz，而对称指标从1（最不对称）变化到5（最对称）。n不必是整数。为了精确拟合的目的，n理论上可以是任何实数。不过在实践中，n通常不应超过5。正如较早前提到的，既然这个专门设计函数的形心频率和对称属性分别由两个不同的参数控制，所以用户利用这个函数精确拟合各种不对称源振幅谱是非常容易的。

根据射线理论，即，将方程式（2）和（3）代入（1）中，如果源振幅谱由加权频率指数函数即S(f)=F(f)近似，那么已接收地震波形的振幅谱可以写成以下形式

$R\left(f\right)=\mathrm{GH}\left(f\right)S\left(f\right)=\mathrm{AG}{f}^n\exp [-f({\∫}_{\mathrm{ray}}\frac{\mathrm{\π}}{\mathrm{Qv}}\mathrm{dl}+\frac{1}{f_0})]---\left(5\right)$

已接收地震波形的形心频率可以被计算为

$f_R=\frac{{\∫}_0^{+\∞}\mathrm{fR}\left(f\right)\mathrm{df}}{{\∫}_0^{+\∞}R\left(f\right)\mathrm{df}}=\frac{n+1}{{\∫}_{\mathrm{ray}}\frac{\mathrm{\π}}{\mathrm{Qv}}\mathrm{dl}+\frac{1}{f_0}}---\left(6\right)$

既然源振幅谱的形心频率f_s是（n+1）f₀，那么源振幅谱与已接收信号振幅谱之间的形心频率迁移可以轻松获得如下

$\mathrm{\Δf}=f_S-f_R=(n+1)(f_0-\frac{1}{{\∫}_{\mathrm{ray}}\frac{\mathrm{\π}}{\mathrm{Qv}}\mathrm{dl}+\frac{1}{f_0}})---\left(7\right)$

沿射线路径的累积衰减现在可以通过求解方程式（7）从形心频率迁移导出：

${\∫}_{\mathrm{ray}}\frac{\mathrm{\π}}{\mathrm{Qv}}\mathrm{dl}=\frac{\mathrm{\Δf}}{f_0[(n+1)f_0-\mathrm{\Δf}]}=\frac{\mathrm{\Δf}}{f_0(f_S-\mathrm{\Δf})}=\frac{\mathrm{\Δf}}{f_0{f}_R}=\frac{(n+1)\mathrm{\Δf}}{f_S{f}_R}---\left(8\right)$

方程式（8）表示利用加权频率指数函数，1/Q通过非常简单的线性关系链接到形心频率迁移。

用于i^th测量（i^th射线）的方程式（8）的离散形式是

$\underset{j}{\mathrm{\Σ}}\frac{\mathrm{\π}{l}_j^i}{Q_j{v}_j}=\frac{\mathrm{\Δ}{f}^i}{f_0{f}_R^i}---\left(9\right)$

其中上角文字i是测量指标以及下角文字j是网格指标，以及是用于第i个测量的在第j个网格中的射线长度。

在收集所有测量后，方程式（9）可以写成矩阵形式

Ax＝d    （10）

在方程式（10）中，A是核心矩阵，其元素由方程式（11）定义

$A_{\mathrm{ij}}=\frac{\mathrm{\π}{l}_j^i}{v_j}---\left(11\right)$

其中x是未知数的矢量，即，

x_j＝1/Q_j    （12）

以及d是由方程式（13）定义的测量矢量

$d_i=\frac{\mathrm{\Δ}{f}^i}{f_0{f}_R^i}---\left(13\right)$

剩下的任务是针对Q求解方程式（10）。既然测量的数据不可避免地被噪音污染，所以方程式（10）的线性***是病态条件化的并且没有唯一的求解办法。因此，这个***可以看成是最小二乘问题，以及用于二次方程式编程问题近似求解办法的一个解决方案

min||Ax-d||    （14）

其中||...||表示欧几里得矢量。

正如较早前所讨论的，为了消除不切实际的负Q值和合并感兴趣领域Q值范围的先验信息，问题（14）可以转换为带有箱约束的最小化问题

min||Ax-d||  条件是1<x<u    （15）

其中l和u是存储Q值的下边界和上边界的矢量。可用于求解优化问题（15）的优选算法类型是例如由Morigi等人（2007）公开的多指标有效集方法。Morigi等人的2007年论文的全部内容在允许其的专利管辖权范围内合并于此，相比于其他约束层析成像算法（Rickett，2006，Delbos等人，2006），这个算法能够更高效并且更有效处理箱约束。这个新型有效集方法的关键特征是其允许一次多指标的有效集的更新。

这个约束优化算法的两级结构化过程可以描述如下。

1.初始化：给出初始模型x⁰和允许误差ε，迭代求解非约束优化问题

||Ax-d||＜ε（16）

以获得非约束解其中||...||表示欧几里得矢量范数，即 $\left|\right|x\left|\right|=\sqrt{({x_1}^2+{x_2}^2+...+{x_n}^2)}.$

2.从k=1开始启动k^th外迭代：将正交投影到可行的矢量集上，以在第一外迭代x^k得到近似的约束解。正交投影由方程式（17）表示

$x_i^k=\left\{\begin{array}{cc}{l}_i& {\hat{x}}_i^k<l_i\\ u_i& {\hat{x}}_i^k>u_i\\ {\hat{x}}_i^k& l_i\&le;{\hat{x}}_i^k\&le;u_i\end{array}\right.---\left(17\right)$

其中上角文字是外迭代指标，以及下角文字是网格指标。满足方程式（17）的首要两个条件中任意一个的这些被称为用于当前迭代的变量的有效集，或者它们的指标被称为指标有效集。也就是说，它们是由箱约束中的一个或另一个影响的变量（指标）。

3.评估余量r^k＝Ax^k-d，如果||r^k||＜ε，那么约束的近似解x^k满足停止标准。终止迭代并且x^k是最终的解。如果不是，则进行步骤4。

4．计算拉格朗日乘子λ^k=A^T(Ax^k–d)并通过设计对角矩阵C^k更新有效集，其中C^k的对角元素定义为

5.通过利用共轭梯度方法，迭代（这是内迭代）求解下列非约束最小化问题，以获得调整矢量y^k。

||B^ky^k+r^k||＜ε（19）

其中B^k＝AC^k。

6.调整解x^k，以通过方程式（20）获得用于下一个外迭代的非约束解

${\hat{x}}^{k+1}=x^k+C^k{y}^k---\left(20\right)$

接着，返回步骤2并进入下一个外迭代。

在本发明的实际应用中，上述方法中六个步骤的某些或全部（通常是全部）利用计算机执行。

上述算法在下列两个方面不同于常规有效集类型方法（Bierlaire等人，1991；1984；Nocedal和Wright，1999）：本发明方法让步骤1和步骤5运行直到满足最优条件（16）和（19）；接着，在步骤4，根据来自以前非约束内循环迭代的是否满足以下两个条件中的一个，通过移除或添加一个或更多（网格）指标，更新有效集的总数：

${\hat{x}}_i^k<l_i$

${\hat{x}}_i^k>u_i$

相反，在常规有效集方法中，一次仅更新有效集中的一个指标。换句话说，近似解一碰到上边界或下边界，常规有效集方法就更新有效集，但是本发明方法直到算法收敛到允许误差才更新有效集。常规有效集算法确保余量是单调递减的；在实际应用中，其示出非常缓慢的收敛速率，这是因为其往往退化到重新启动的最速下降法，尤其是在矩阵是高度病态时。这个新类型的有效集算法不保证给出单调递减数据失配。不过，经验表明，如果选择恰当允许误差，则这个算法在快速收敛速率方面完成的非常好并且产生满足允许误差的可行解。

如同上述步骤1-6所例示的，本发明有效集类型算法的另一个优势是用于启动优化的开始的Q模型可以是任意的，这意味着开始的Q模型可以在上或下边界或甚至在边界之外设置。这个特征在某些情况下是非常有用的。

本发明可以根据图3所示流程图实施。如同在常规的基于射线的层析成像算法中一样，在步骤60，速度模型20和源/接收器位置30被输入到射线追迹代码中，以输出射线路径信息，其用于在步骤70中构建核心矩阵（在方程式10和11中的A）。在步骤40，地震数据被分析并且源振幅谱由加权频率指数函数（方程式3）近似，在此期间对称指标n和特征频率f₀被确定。接着在步骤50计算所有地震轨迹的形心频率，以获得形心频率迁移（方程式7）。在步骤80和步骤90，收集先验信息以便建造开始的Q模型和Q范围地图。其后，在步骤100，所有地震轨迹、核心矩阵、开始的Q模型以及Q范围地图的形心频率迁移被输入到多指标有效集约束优化代码中（上述6步骤算法），以执行Q模型重构。在步骤110，重构的Q模型由用户判定。如果重构的Q模型被接受，则该Q层析成像过程结束。否则，用户重新建造开始的Q模型和/或Q范围地图，并再次实施约束的优化，直到重构的Q模型是令人满意的。上述的括号***内容指的是本发明的示例性实施例。本文的描述集中在步骤400和步骤100，其是本发明的主要要点。除非另有说明，在图3中示出的其他步骤在Q层析成像领域是众所周知的，而且为了本发明的目的，其可以以任何已知或后来开发的方式执行。而且，虽然既实施本发明对步骤40的改善又实施本发明对步骤100的改善是优选的，但是可以在没有实施另一个的情况下实施其中一个，即标准的途径可用于履行步骤40或步骤100的功能。

可以添加各种可选步骤或改进到图3所示的基本步骤。例如：

1）所有地震轨迹可以在步骤40之前被预处理：第一批到达者被分窗口并被隔离，计算隔离子波的振幅谱，以及在振幅谱上实施多项式拟合。（原始地震数据往往是噪声非常多的，因此，不是利用高斯函数尝试拟合原始地震振幅谱，在实践中的第一步可以是原始地震振幅谱的多项式拟合，以将其转换为平滑曲线。在应用本发明方法中，这个相同中间步骤可能是有用的，这之后，用加权频率指数函数拟合平滑多项式曲线就更加容易。）

2）可以添加正则矩阵，以稳定步骤100的优化过程。

3）可以从一个或更多速度异常建造开始的Q模型（步骤80）。

4）当地震轨迹数量太小时，可以采用数据规则化技术。

5）在步骤100期间，可以利用共轭梯度方法或LSQR方法（共轭梯度方法的变种；参见Paige和Saunders，1982）实施约束优化算法的内迭代。

6）在步骤40期间，如果源谱不可获得，那么可以从接收的地震波形估算它，接收的地震波形的关联射线不穿越带有可能的Q异常的区域。

示例

在这部分，呈现Q层析成像的合成示例。本发明的Q层析成像算法的性能与利用带有或不带有箱约束的高斯函数的常规形心频率迁移Q层析成像算法的性能作比较。图4示出不对称的源振幅谱。源振幅谱的形心频率可以被数值地计算为55.7Hz，而峰值频率是33Hz。通过高斯函数源振幅谱拟合，高斯函数的峰值频率是55.7Hz而方差是1246Hz²。通过加权频率指数函数拟合，在方程式（3）中的参数n和f₀分别是1.45和22.7Hz。图5A是速度模型和通过运行射线追迹代码获得的射线路径，而图5B是真实的1/Q模型。在这个合成数据测试中，使用了50个源。对于每个源，存在90个接收器。因此，总共有4500个轨迹。图5A-B到8A-B是其中绘出量的显示黑白复制品，根据绘图中的颜色标度，用颜色表示的深度和距离。

图6A是利用带有高斯函数拟合和使用非约束优化的常规形心频率迁移Q层析成像方法的重构的1/Q模型。因此，由于由高斯源振幅谱拟合引入的大量误差，重构的1/Q模型明显偏离真实的1/Q模型。图6B示出重构的1/Q模型与真实1/Q模型之间的差异。接着，高斯函数由本发明的加权频率指数函数替换，以拟合源振幅谱和重构的1/Q模型，其中图7A中示出非约束优化算法。相比于利用高斯拟合获得的结果，结果极大改善了。不过，在重构的1/Q模型中存在许多人工产物，并且在某些区域中的Q值是负的，这在物理上是不真实的。图8A示出利用本发明的加权频率指数函数拟合的Q层析成像结果，这次也使用本发明带有箱约束的优化算法，其强迫重构的1/Q值在0到0.05之间的范围内。图8B示出重构的1/Q模型与真实1/Q模型之间的差异。这个结果比之前两个更好，这是因为在该结果中没有负的Q值以及存在很少人工产物。

在本发明的某些实施例中，通过计算不穿越示出Q异常的地下区域的地震轨迹的平均振幅谱，执行对源振幅谱的估算。平滑可以应用于平均振幅谱，其可以利用多项式拟合技术完成。

在本发明的某些实施例中，可以从地震数据轨迹的第一批到达者计算接收信号的振幅谱，这些地震数据轨迹可以通过窗口被隔离。

在本发明的某些实施例中，可以从地下声学速度模型建造开始的1/Q模型。可以通过线性映射从地下声学速度模型中的一个或更多异常建造开始的1/Q模型。

在本发明的某些实施例中，1/Q范围地图可以具有都大于0的下边界。

前面的描述以本发明特定实施例为例子进行说明，不过对于本领域的技术人员显而易见的是，针对本文所述实施例的许多修改和变化是可能的。正如所附权利要求定义的一样，所有这样的修改和变化都应当在本发明的范围内。

参考文献

Bierlaire,M.,Ph.L.Toint,和D.Tuyttens,“On iterative algorithms forlinear least squares problems with bound constraints,”Linear AlgebraAppl 143,111-143(1991).

Delbos,D.,J.Ch.Gilbert,R.Glowinski,和D.Sinoquet,“Constrainedoptimization in seismic reflection tomography:a Gauss-Newtonaugmented Lagrangian approach,”Geophysics 164,670-684(2006).

Gao,F.,G.Fradelizio,A.Levander,G.Pratt,和C.Zelt,“Seismic velocity,Q,geological structure and lithology estimation at a ground watercontamination site,”75th SEG，Expanded Abstracts,1561-1564,Soc.ofExpl.Geophys(2005).

Hicks,G.J.和R.G.Pratt,“Reflection waveform inversion using localdescent methods:Estimating attenuation and velocity over a gas-sanddeposit，”Geophysics 66,598–612(2001).

Liao,Q.和G.A.McMechan,“Multifrequency viscoacoustic modelingand inversion,”Geophysics 61,No.5,1371-1378(1996).

P.,“Solving the minimal least squares problem subject to boundson the variables,”BIT 24(1984).

Morigi,S.,L.Reichel,F.Sgallari,和F.Zama,“An iterative method forlinear discrete ill-posed problems with box constraints,”Journal ofComputational and Applied Mathematics 198,505-520(2007).

Nocedal,J.和S.J.Wright,Numerical Optimization,Springer,New York(1999).

Paige,C.C.和M.A.Saunders,“LSQR:An algorithm for sparse linearequations and sparse least squares,”ACM Trans.Math.Software 8,43-71(1982).

Plessix,R.-E.,“Estimation of velocity and attenuation coefficient mapsfrom crosswell seismic data,”Geophysics 71,S235–S240(2006).

Pratt,R.G.,K.Bauer，和M.Weber,“Crosshole waveform tomographyvelocity and attenuation images of arctic gas hydrates,”73rd  SEG，Expanded Abstracts,Society of Exploration Geophysics,2255-2258(2003).

Rickett,J.E.,“Method for estimation of interval seismic quality factor,”U.S.Patent No.7,376,517(2006).

Rossi,G.,D.Gei,G.G.Madrussani,和J.M.Carcione,“Attenuation tomography:An application to gas-hydrate and free-gasdetection,”Geophysics 55,655-669(2007).

Spencer，T.W.,J.R.Sonnad,和T.M.Butler，“Seismic Q-stratigraphy ordissipation,”Geophysics 47,16-24(1982).

Tonn,R.,“The determination of the seismic quality factor Q from VSP data:A comparison of different computational method,”GeophysicalProspecting 39,1–27(1991).

Quan,Y.和J.M.Harris,“Seismic attenuation tomography using thefrequency shift method,”Geophysics 62,895-905(1997).

Watanabe,T.,K.T.Nihei,S.Nakagawa,和L.R.Myer，“Viscoacousticwaveform inversion of transmission data for velocity and attenuation,”J.Acoust.Soc.Am.115,3059–3067(2004).

Wright,S.J.,“Primal-Dual Interior Point Methods,”SIAM(1997).

Wu,H.,和J.M.Lees,“Attenuation structure of Coso geothermal data,California from wave pulse widths,”Bulletin of the Seismological Societyof America 86,1574-1590(1996).

Zhang,C.,“Seismic Absorption Estimation and Compensation,”PhD thesis,University of British Columbia(2008).

Claims (20)

1.一种基于射线的形心频率迁移Q层析成像方法，用于从接收器在利用地震源的勘探中测量的地震数据重构地下Q模型，所述方法包括选择数学函数来近似所述地震源的振幅谱，以便计算由于地面衰减而造成的所述谱的形心频率迁移，并使所述形心频率迁移与由品质因数Q的倒数表示的衰减相关，以及利用计算机通过迭代、线性优化求解Q或1/Q，其中所述优化具有箱约束，以将估算的Q值保持在由上边界和下边界规定的与位置有关的范围内。

2.根据权利要求1所述的方法，其中所述具有箱约束的优化通过多指标有效集方法求解，所述多指标有效集方法允许一次多网格指标的有效集更新，其中网格指标表示地下位置。

3.根据权利要求1所述的方法，其中所选择的数学函数是频率的加权频率指数函数。

4.根据权利要求3所述的方法，其中所述频率的加权频率指数函数具有两个参数，所述两个参数被调整以提供到所述地震源的振幅谱的拟合。

5.根据权利要求4所述的方法，其中所述两个参数是用于带宽控制的特征频率和对称指标，每个是正实数。

6.根据权利要求5所述的方法，其中所述频率的加权频率指数函数能够表述成下面的形式

$F\left(f\right)={\mathrm{Af}}^n\exp (-\frac{f}{f_0})$

其中f是频率，A是用于振幅换算的常数，以及f₀是所述特征频率，以及n是所述对称指标。

7.根据权利要求1所述的方法，其进一步包括：

(a)估算所述源的振幅谱和计算其形心频率；

(b)通过频率的加权频率指数函数来近似所述源的振幅谱；

(c)计算所述地震数据的轨迹的第一批到达者的振幅谱；

(d)计算形心频率迁移，其是在(c)中计算的振幅谱的形心频率与所述源的振幅谱的计算形心频率之间的差异；

(e)以所述形心频率迁移和在(c)中计算的振幅谱的形心频率建造测量矢量d；

(f)利用地下声学速度模型和来自所述勘探的源接收器信息，在计算机上运行射线追迹代码；

(g)以射线段长度和相应声学速度建造核心矩阵A；

(h)从可用信息建造地下的开始Q模型，所述开始模型为所述开始模型中的每个单元规定1/Q的值；

(i)生成贯穿所述模型为1/Q提供箱约束的1/Q范围地图，所述箱约束基于可用信息；

(j)执行迭代优化，其中计算机为矢量x的分量求解问题min||Ax-d||，服从箱约束，其中x_j＝1/Q_j，指标j表示所述模型中

的第j个单元，从而重构1/Q值的体积作为地下的深度和横向位置的

函数。

8.根据权利要求7所述的方法，其中服从箱约束的所述迭代优化通过采用一类型的有效集方法执行，所述有效集方法在所述优化期间更新有效集的多个指标。

9.根据权利要求8所述的方法，其中所述迭代优化具有外迭代循环和内迭代循环，以及所述内迭代循环执行确定对x_i的调整的非约束优化，以使||Ax-d||最小。

10.根据权利要求9所述的方法，其中在每个内循环迭代后产生针对每个模型单元i的x_i，通过以下步骤开始下一个外循环：对照箱约束测试x_i和约束违反所述约束的那些x_i，所述受约束的x_i称为所述有效集，接着测试||Ax-d||是否小于ε，如果不是，则进行下一个内迭代。

11.根据权利要求10所述的方法，其中直到所述内循环的非约束迭代收敛以满足选择的最优条件，才更新所述有效集的总体。

12.根据权利要求2所述的方法，其中所述多指标有效集方法使用共轭梯度求解程序或LSQR求解程序。

13.根据权利要求6所述的方法，其中所述使所述形心频率迁移Δf与1/Q相关可以数学表示为

$\mathrm{\&Delta;f}=f_S-f_R=(n+1)(f_0-\frac{1}{{\&Integral;}_{\mathrm{ray}}\frac{\mathrm{\&pi;}}{\mathrm{Qv}}\mathrm{dl}+\frac{1}{f_0}})$

其中f_S和f_R分别是所述地震源的振幅谱和接收器检测的振幅谱的形心频率，v是声学速度，以及dl是射线路径的增量。

14.根据权利要求1所述的方法，其进一步包括使用在地震成像中的Q或1/Q的求解值，用于碳氢化合物勘探或用于碳氢化合物储集层特性描述。

15.一种基于射线的形心频率迁移Q层析成像方法，用于从接收器在利用地震源的勘探中测量的地震数据重构地下Q或1/Q模型，所述方法包括利用频率的加权频率指数函数来近似所述地震源的振幅谱，以便计算由于地面衰减而造成的所述谱的形心频率迁移，并使所述形心频率迁移与由品质因数Q的倒数表示的衰减相关，以及通过利用计算机执行的迭代、线性优化求解Q或1/Q。

16.根据权利要求15所述的方法，其中所述优化具有箱约束，以将估算的Q值保持在由上边界和下边界规定的与位置有关的范围内。

17.根据权利要求15所述的方法，其中所述具有箱约束的优化通过多指标有效集方法求解，所述多指标有效集方法允许一次多网格指标的有效集更新，其中网格指标表示地下位置。

18.根据权利要求15所述的方法，其中所述频率的加权频率指数函数具有两个参数，所述两个参数被调整以提供到所述地震源的振幅谱的拟合，其中所述两个参数是用于带宽控制的特征频率和对称指标，每个是正实数。

19.根据权利要求18所述的方法，其中所述频率的加权频率指数函数可以表述成下面的形式

$F\left(f\right)={\mathrm{Af}}^n\exp (-\frac{f}{f_0})$

其中f是频率，A是用于振幅换算的常数，f₀是所述特征频率，以及n是所述对称指标。

20.根据权利要求15所述的方法，其进一步包括使用在地震成像中的Q或1/Q的求解值，用于碳氢化合物勘探或用于碳氢化合物储集层特性描述。