Q层析成像方法
阅读：457发布：2020-05-11
专利汇可以提供Q层析成像方法专利检索，专利查询，专利分析的服务。并且本 发明 涉及一种通过执行基于射线（60）的形心 频率 迁移（50）Q 层析成像 ，从 地震 数据（10）重构地下Q模型（110）的方法。地 震源 波形 的振幅谱由频率（40）的加权频率指数函数近似，具有两个参数以便调整以拟合频率迁移数据，从而提供对各种不对称源振幅谱的拟合。箱约束可用于优化程序，以及用于速度层析成像的多指标有效集方法是用于实施箱约束（100）的优选技术。，下面是Q层析成像方法专利的具体信息内容。
1.一种基于射线的形心频率迁移Q层析成像方法，用于从接收器在利用地震源的勘探中测量的地震数据重构地下Q模型，所述方法包括选择数学函数近似所述地震源的振幅谱，以便计算由于地面衰减而造成的所述谱的形心频率迁移，并使所述形心频率迁移与由品质因数Q的倒数表示的衰减相关，以及利用计算机通过迭代、线性优化求解Q或1/Q，其中所述优化具有箱约束，以将估算的Q值保持在由上边界和下边界规定的与深度有关范围内的。
2.根据权利要求1所述的方法，其中所述具有箱约束的优化通过多指标有效集方法求解，所述多指标有效集方法允许一次多网格指标的有效集更新，其中网格指标表示地下位置。
3.根据权利要求1所述的方法，其中所选择的数学函数是频率的加权频率指数函数。
4.根据权利要求3所述的方法，其中所述频率的加权频率指数函数具有两个参数，所述参数被调整以提供到所述地震源的振幅谱的拟合。
5.根据权利要求4所述的方法，其中所述两个参数是用于带宽控制的特征频率和对称指标，每个是正实数。
6.根据权利要求5所述的方法，其中所述频率的加权频率指数函数可以表述成下面的形式
其中f是频率，A是用于振幅换算的常数，以及f0是所述特征频率，以及n是所述对称指标。
7.根据权利要求1所述的方法，其进一步包括：
（a）估算所述源的振幅谱和计算其形心频率；
（b）通过频率的加权频率指数函数近似所述源的振幅谱；
（c）计算所述地震数据的轨迹的第一批到达者的振幅谱；
（d）计算形心频率迁移，其是在（c）中计算的振幅谱的形心频率与所述源的振幅谱的计算形心频率之间的差异；
（e）以所述形心频率迁移和在（c）中计算的振幅谱的形心频率建造测量矢量d；
（e）利用地下声学速度模型和来自所述勘探的源接收器信息，在计算机上运行射线追迹代码；
（f）以射线段长度和相应声学速度建造核心矩阵A；
（g）从可用信息建造地下的开始Q模型，所述开始模型为所述开始模型中的每个单元规定1/Q的值；
（h）生成贯穿所述模型为1/Q提供箱约束的1/Q范围地图，所述箱约束基于可用信息；
（i）执行迭代优化，其中计算机为矢量x的元素求解问题min||Ax-d||，服从箱约束，其中xj＝1/Qj，指标j表示所述模型中的第j个单元，从而重构1/Q值的体积作为地下的深度和横向位置的函数。
8.根据权利要求7所述的方法，其中服从箱约束的所述迭代优化通过采用一类型的有效集方法执行，所述有效集方法在所述优化期间更新有效集的多个指标。
9.根据权利要求8所述的方法，其中所述迭代优化具有外迭代循环和内迭代循环，以及所述内迭代循环执行确定对xi的调整的非约束优化，以使||Ax-d||最小。
10.根据权利要求9所述的方法，其中在每个内循环迭代后产生针对每个模型单元i的xi，通过以下步骤开始下一个外循环：对照箱约束测试xi和约束违反所述约束的那些xi，所述受约束的xi称为所述有效集，接着测试||Ax-d||是否小于ε，如果不是，则进行下一个内迭代。
11.根据权利要求10所述的方法，其中直到所述内循环的非约束迭代收敛以满足选择的最优条件，才更新所述有效集的总体。
12.根据权利要求2所述的方法，其中所述多指标有效集方法使用共轭梯度求解程序或LSQR求解程序。
13.根据权利要求6所述的方法，其中所述使所述形心频率迁移Δf与1/Q相关可以数学表示为
其中fS和fR分别是所述地震源的振幅谱和接收器检测的振幅谱的形心频率，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是频率，A是用于振幅换算的常数，f0是所述特征频率，以及n是所述对称指标。
20.根据权利要求15所述的方法，其进一步包括使用在地震成像中的Q或1/Q的求解值，用于碳氢化合物勘探或用于碳氢化合物储集层特性描述。
Q层析成像方法
[0001] 相关申请的交叉引用
技术领域
背景技术
[0004] 针对岩石属性的特性描述和适当的振幅随偏移距变化（AVO）分析，需要考虑地震衰减效应。在迁移过程中，需要地震衰减信息来补偿吸收效应，以提高迁移图像质量。因此，地震衰减的估算对于储油层的检测和监视是必要的。
[0005] 地震衰减可以由品质因数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层析成像类别。
Hicks和Pratt，2001；Pratt等人，2003；Watanabe等人，2004；Gao等人，2005）。基于波方程的层析成像在物理上更为精确但是计算昂贵并且对3D情况不适用。本发明属于基于射线的Q层析成像类别。
[0006] Q层析成像的一个主要问题是如何以最小的近似值和最大的灵活性建立Q模型与地震数据之间的联系。广泛使用的方法是基于Q与地震振幅衰减之间的关系。另一个方法使用地震形心频率降档来估算品质因数Q。后面的这个方法被认为更加稳健，这是因为这个方法与几何扩散效应和反射/传输损失无关。不过，常规形心频率迁移方法仅能使用高斯函数、矩形波函数或三角形函数来拟合源振幅谱，这会引入明显的误差，这是因为在多数情况下，源谱不能由这些函数近似。本发明包括加权频率指数函数，其被设计为拟合各种不对称源振幅谱，从而通过大大减少源振幅谱拟合误差，改善Q层析成像的精度。
[0007] 在多数现有Q层析成像算法中，优化部分是基于无约束的优化方法或基于简单的非负约束的优化方法。因此，这些Q层析成像算法花费大量计算时间或产生许多人工产物和不合实际的Q模型（例如，负的Q值或极端低的Q值），尤其是当地震数据被噪音污染时。包括带有箱约束的高效优化算法的本发明能够改善重构Q模型的质量和可靠性。下面进行现有技术的更详细讨论。
[0008] 地震衰减层析成像（Q层析成像）已被研究许多年并取得很大进展。基于射线Q层析成像算法的两个主要组成部分是1）用于构建Q层析成像的数学模型的简单但是精确的地震数据与Q模型之间的关系；2）解这个数学问题的可靠且稳健的优化算法。已经开发或提议用于构建这两个组成部分的许多技术。这些技术讨论如下。
[0009] 建立地震数据与Q模型之间的联系
[0010] 估算Q的最简单直接的方法是谱比方法（Spencer等人，1982；Tonn，1991），其中两个地震波形之间谱比的对数被计算为频率的函数，并且这个函数由线性频率函数近似，其斜率被视为累积的地震衰减并最终与沿波传播路径的Q值相关。理想地，这个方法假设几何扩散效应和反射/传输损失效应与频率无关，消除了几何扩散效应和反射/传输损失效应。在实际应用中，由于子波重叠、线性拟合的不确定性以及许多其他因素，这个方法是相对不可靠的。
[0011] Rickett（2006）提出利用时间-频率分析技术辅助的谱比层析成像扩展方法。这个方法被宣称对绝对比例不敏感并应用于利用垂直地震剖面（VSP）的Q剖面估算的应用中。在这个方法中，描述与频率无关的振幅变化的对数-振幅标量包含在未知数中，这大大增加了未知数的数量并降低算法的效率。而且，在利用表面地震波反射数据的2D/3DQ层析成像中，对数-振幅标量不仅是位置的函数，而且是射线的函数，这使过程非常复杂。
[0012] 基于地震子波上升时间的变化与沿传播路径的1/Q剖面是线性相关的，Wu和Lees（1996）报告了利用地震学中的上升时间的地震衰减层析成像方法。遗憾的是，这个方法在勘探地球物理中是不实用的，这是因为子波不可避免地被噪音、散射效应、重叠等污染。
[0013] 应当指出，地震子波振幅谱的形状几乎完全受品质因数Q影响，并且开发了用于Q估算的峰值频率变化方法（Zhang，2008）。这个方法是有吸引力的，但是在实际中，存在精确测量峰值频率变化的困难。而且，由于仅使用在单个频率的信息，所以Q估算的不确定性大。
[0014] 更稳健的方法由Quan和Harris引入（1997），其中，地震波形的整个频段的信息被用于计算形心频率降档，接着通过简单封闭形式公式，使形心频率迁移与沿射线路径的Q剖面相关。这个方法本质上是完全不受几何扩散和反射/传输损失影响的。这个方法的局限在于源振幅谱必须是高斯、矩形波或三角形函数。众所周知的是地震振幅谱从来都不是矩形波或三角形函数。而且，其通常是不对称的并且与高斯函数有很大不同。若这个不对称振幅谱由高斯函数近似，则会在Q模型的重构中引入明显误差。因此，若在不丢失记录的地震数据的形心频率与沿射线路径的Q剖面之间关系的简单性质的情况下，存在能够用于更精确拟合各种地震频谱的函数，则这个稳健方法在实际Q层析成像应用中会更加精确。
[0015] 用于Q层析成像的约束优化算法
[0016] 当建立地震数据与Q模型之间的关系时，基于射线的Q层析成像问题可以由线性优化问题描述。在大多数现有Q层析成像算法中，这个线性优化问题通过利用Krylov子空间方法迭代求解，例如共轭梯度方法和LSQR方法，而不应用任何约束（Quan和Harris，1997；Plessix2006，Rossi等人，2007）。只要地震数据具有高信噪比（SNR），那么这些算法工作得很好。不过，地震数据从来都不是干净的；在处理真实的场数据时，这些非约束的优化算法总是产生某些负的Q值或极端小的正Q值，这些是物理上不现实的。而且，在某些情况下，Q值范围的先验信息是已知的。在这些情况下，这个信息需要通过箱约束包含在层析成像算法中，以便提供更加可靠的Q层析成像结果。
[0017] Rickett（2006）开发出具有约束的Q估算算法。但是他的算法采用非负约束，而不是箱约束，这意味着消除了负的Q值，但仍然存在极低的正Q值。Rickett（2006）报告了应用非负约束的两种方法。第一种方法是非线性变换方法。在这个非线性变换方法中，y变量Q由e 代替，并且求解y而不是1/Q。通过这样做，Q被迫为正，但是整个系统会是非常非线性的。为了实现这个目标，由此产生的优化系统利用高斯-牛顿途径求解，而这是很昂贵的。这个方法的另一个缺点是，在优化过程中当Q值很大时，基于梯度的优化算法将是停滞的或收敛非常缓慢，即Q值会保持在那里并且不再变化。在最恶劣情况下，若在开始模型中的y值是无限的，那么成本函数的梯度是0并且优化算法不执行搜索。应用Rickett（2006）报告的非负约束的第二个方法是通过平滑技术，强制衰减的单调递增特性。这个方法利用VSP数据对Q估算是有效起作用的，但在利用表面反射地震数据的2D Q层析成像中会失效。
[0018] 本发明人知道，不存在带有箱约束的现有Q层析成像算法来实施在上边界和下边界规定范围内的估算Q值。不过，在某些其他地球物理应用中，例如速度层析成像中采用带有箱约束的优化算法。例如，Delbos等人（2006）开发出具有箱约束的地震反射层析成像算法。在他们的算法中，利用高斯—牛顿增强朗格拉日方法求解约束优化问题，而关联的朗格拉日问题——另一个约束优化问题，通过梯度投影方法、有效集方法以及共轭梯度方法的组合来求解。他们使用的有效集方法是常规的，并且由于该算法更新有效集，一次一个约束，因此是低效的（Bierlaire等人，1991； 1984；Nocedal和Wright，1999）。当箱约束的数量巨大时，该算法的收敛速率可以非常慢。
[0019] 在本发明中，在数学领域中的最近进展（其被称为多指标有效集方法（Morigi等人，2007））被用于执行带有箱约束的Q层析成像，相比于非约束的Q层析成像算法和利用常规有效集方法的约束算法，其明显改善在Q重构质量和算法效率方面的Q层析成像算法的性能。
发明内容
[0020] 在一个实施例中，本发明是基于射线的形心频率迁移Q层析成像方法，其从在利用地震源的勘测中由接收器测量的地震数据重构1/Q的地下深度模型，其包括选择数学函数来近似地震源的振幅谱，以便计算由于地球衰减引起的谱的形心频率迁移，并且将所述形心频率迁移与由品质因数Q的倒数表示的衰减相关，以及通过迭代的线性优化求解Q或1/Q，其中所述优化具有箱约束，以将估算的Q值保持在由上边界和下边界规定的与位置相关的范围内。该约束优化可以由多指标有效集方法求解，该方法允许一次多网格指标的有效集更新，其中网格指标表示地下位置。选择的数学函数可以是频率的加权频率指数函数。
无论在优化中是否使用箱约束，使用这个函数来近似地震源的振幅谱将是有利的。
无论在优化中是否使用箱约束，使用这个函数来近似地震源的振幅谱将是有利的。
[0021] 除了其他优点之外，由本发明方法产生的重构的地下Q模型可以有利地用于地震成像中，以补偿振幅模糊、频率损失以及由粘声波覆盖层，例如天然气储集层导致的相位失真效应。通过在地震成像过程中包含由本发明提供的更精确Q模型，能够明显改善地质结构图像的质量。此外，由于Q对某些岩石和流体属性，例如流体饱和度与孔隙度很敏感，所以重构的Q模型在储集层特征描述的应用中是有益的。在这样的应用中，本发明成为碳氢化合物勘探、开发或生产的方法。
[0023] 通过参考下列详细描述及其随附绘图，可以更好理解本发明及其优势，其中：
[0024] 图1示出本发明至少某些实施例中设计用于拟合各种不对称源振幅谱的加权频率指数函数的带宽变化，针对特性频率参数的各种值和称为对称指标的第二参数的选择值在n=2时的情况；
[0025] 图2示出加权频率指数函数 在n=1，2，3，4，5而形心频率（n+1）f0固定在30Hz时的曲线图；
[0026] 图3示出在Q层析成像中实施本发明的基本步骤的流程图；
[0027] 图4到8A-B属于合成示例应用，其中
[0028] 图4示出源振幅谱；
[0029] 图5A-B示出速度模型和射线路径（5A）以及真实的1/Q模型（5B）；
[0030] 图6A示出利用具有高斯函数拟合和没有箱约束的常规形心频率迁移Q层析成像方法的重构的1/Q模型；以及图6B示出图6A的重构1/Q模型与图5B的真实1/Q模型之间的差异；
[0031] 图7A示出利用带有本发明的加权频率指数函数拟合、但是没有箱约束的形心频率迁移Q层析成像方法的重构的1/Q模型；图7B示出图7A的重构1/Q模型与图5B的真实1/Q模型之间的差异；以及
[0032] 图8A示出利用带有本发明的加权频率指数函数拟合和根据本发明的箱约束的形心频率迁移Q层析成像方法的重构的1/Q模型；图8B示出图8A的重构1/Q模型与图5B的真实1/Q模型之间的差异。
[0034] 将结合示例性实施例来描述本发明。不过在某种程度上，下列详细描述具体到本发明的特定实施例或特定用途，这仅仅是为了说明目的，而不应解释为对本发明范围的限制。相反，正如所附权利要求所定义的，意欲覆盖包含在本发明范围内的所有替换、修改和等同物。
具体实施方式
[0035] 本发明包括从地震数据重构2D或3D地震品质因数（Q）模型的方法，称为Q层析成像的技术领域。
[0036] 本发明在至少某些实施例中的主要特征如下。专门设计的函数，加权频率指数函数分析和拟合源子波振幅谱。通过调整加权频率指数函数的两个参数实施源子波振幅谱拟合。已接收地震波形的形心频率相对源子波形心频率的迁移被计算并输入到带有箱约束的优化算法中，以便重构Q模型，其中Q值的范围由先验信息预先确定。不像最初的形心频率迁移方法（Quan和Harris，1997），利用专门设计的振幅谱拟合函数，这个Q层析成像算法能够处理非高斯且非对称源振幅谱，以减少由源振幅谱拟合中的不匹配引入的误差。利用这个专门设计的源谱拟合函数，Q层析成像问题转换为带有箱约束的约束优化问题。这个约束优化通过采用多指标有效集方法（Morigi等人，2007）求解，其进一步改善这个Q层析成像算法的精度和稳健性，而不牺牲高效率特征。术语有效集指的是指示的不能够在迭代周期结束时更新的正被优化未知量组的子集，这是因为他们碰到约束，无论是上限还是下限。
[0037] 下一步解释本发明的某些基本理论。
[0038] 如果源子波的振幅谱是S(f)，那么接收的地震波形的振幅谱R(f)可以表述为（Quan和Harris，1997）
[0039] R(f)＝GH(f)S(f) (1)
[0040] 其中G是包括几何扩散效应、反射/传输系数等的、与频率无关的因数。H(f)是描述地震衰减效应的脉冲响应函数，其被公式化为
[0041]
[0042] 其中Q是品质因数以及v是地震（即，声学）波速度。
[0043] 在用于Q估算的形心频率迁移方法中，关键部分是使用分析函数来拟合源振幅谱S(f)，接着导出已接收地震的形心频率与沿波传播路径的Q剖面之间的显性关系。在最初的形心频率迁移方法中，这个显性关系仅通过假设源振幅谱是高斯函数或矩形波函数（除了单个间隔处是等于常数以外，在整个实线上是零的函数）或三角形函数导出。这是传统形心频率迁移方法的主要不足，其会导致Q估算中的大量误差（Rickett，2006；Zhang，2008）。
[0044] 本发明至少某些实施例的一个主要特征是加权频率指数函数被设计用于更精确拟合各种不对称源振幅谱，而不丢失形心频率迁移与沿射线路径的Q剖面之间的简单封闭形式关系。而且，利用这个专门设计的函数实施源幅度谱拟合是非常方便的，这是因为函数的形状和带宽由两个参数分别确定。加权频率指数函数的公式是
[0045]
[0046] 其中A是振幅比例的常数，f0称为特征函数，以及n称为对称指标。特征频率f0是用于带宽控制的参数。如果对称指标n固定，则这个函数的带宽随着特征频率增加而扩展。实际上，F(f)的形心频率具有非常简单的形式，给出如下：
[0047]
[0048] 图1示出当的n=2时F(f)的带宽变化，具有范围从10Hz到50Hz的增量为10Hz的五个不同特征频率。在图1中，当f0从50Hz变化到10Hz时，F(f)的形状保持相对不变，而带宽如预计的那样萎缩。另一方面，对称指标n被用于控制F(f)的对称属性。对称指标n越大，函数F(f)的形状越对称，如图2所示，其中形心频率固定为30Hz，而对称指标从1（最不对称）变化到5（最对称）。n不必是整数。为了精确拟合的目的，n理论上可以是任何实数。不过在实践中，n通常不应超过5。正如较早前提到的，既然这个专门设计函数的形心频率和对称属性分别由两个不同的参数控制，所以用户利用这个函数精确拟合各种不对称源振幅谱是非常容易的。
[0049] 根据射线理论，即，将方程式（2）和（3）代入（1）中，如果源振幅谱由加权频率指数函数即S(f)=F(f)近似，那么已接收地震波形的振幅谱可以写成以下形式
[0050]
[0051] 已接收地震波形的形心频率可以被计算为
[0052]
[0054]
[0055] 沿射线路径 的累积衰减现在可以通过求解方程式（7）从形心频率迁移导出：
[0056]
[0057] 方程式（8）表示利用加权频率指数函数，1/Q通过非常简单的线性关系链接到形心频率迁移。th th
[0058] 用于i 测量（i 射线）的方程式（8）的离散形式是
[0059]
[0060] 其中上角文字i是测量指标以及下角文字j是网格指标，以及 是用于第i个测量的在第j个网格中的射线长度。
[0061] 在收集所有测量后，方程式（9）可以写成矩阵形式
[0062] Ax＝d （10）
[0063] 在方程式（10）中，A是核心矩阵，其元素由方程式（11）定义
[0064]
[0065] 其中x是未知数的矢量，即，
[0066] xj＝1/Qj （12）
[0067] 以及d是由方程式（13）定义的测量矢量
[0068]
[0069] 剩下的任务是针对Q求解方程式（10）。既然测量的数据不可避免地被噪音污染，所以方程式（10）的线性系统是病态条件化的并且没有唯一的求解办法。因此，这个系统可以看成是最小二乘问题，以及用于二次方程式编程问题近似求解办法的一个解决方案[0070] min||Ax-d|| （14）
[0071] 其中||...||表示欧几里得矢量。
[0072] 正如较早前所讨论的，为了消除不切实际的负Q值和合并感兴趣领域Q值范围的先验信息，问题（14）可以转换为带有箱约束的最小化问题
[0073] min||Ax-d|| 条件是1
[0074] 其中l和u是存储Q值的下边界和上边界的矢量。可用于求解优化问题（15）的优选算法类型是例如由Morigi等人（2007）公开的多指标有效集方法。Morigi等人的2007年论文的全部内容在允许其的专利管辖权范围内合并于此，相比于其他约束层析成像算法（Rickett，2006，Delbos等人，2006），这个算法能够更高效并且更有效处理箱约束。这个新型有效集方法的关键特征是其允许一次多指标的有效集的更新。
[0075] 这个约束优化算法的两级结构化过程可以描述如下。
[0076] 1.初始化：给出初始模型x0和允许误差ε，迭代求解非约束优化问题[0077] ||Ax-d||＜ε（16）
[0078] 以 获 得 非 约 束 解 其 中||...|| 表 示 欧 几 里 得 矢 量 范 数，即[0079] 2.从k=1开始启动kth外迭代：将 正交投影到可行的矢量集上，以在第一外迭代kx 得到近似的约束解。正交投影由方程式（17）表示
[0080]
[0081] 其中上角文字是外迭代指标，以及下角文字是网格指标。满足方程式（17）的首要两个条件中任意一个的这些 被称为用于当前迭代的变量的有效集，或者它们的指标被称为指标有效集。也就是说，它们是由箱约束中的一个或另一个影响的变量（指标）。
[0082] 3.评估余量rk＝Axk-d，如果||rk||＜ε，那么约束的近似解xk满足停止标准。终止迭代并且xk是最终的解。如果不是，则进行步骤4。
[0083] 4．计算拉格朗日乘子λk=AT(Axk–d)并通过设计对角矩阵Ck更新有效集，其中Ck的对角元素 定义为
[0084]
[0085] 5.通过利用共轭梯度方法，迭代（这是内迭代）求解下列非约束最小化问题，以获k得调整矢量y。
[0086] ||Bkyk+rk||＜ε（19）
[0087] 其中Bk＝ACk。
[0088] 6.调整解xk，以通过方程式（20）获得用于下一个外迭代的非约束解[0089]
[0090] 接着，返回步骤2并进入下一个外迭代。
[0091] 在本发明的实际应用中，上述方法中六个步骤的某些或全部（通常是全部）利用计算机执行。
[0092] 上述算法在下列两个方面不同于常规有效集类型方法（Bierlaire等人，1991；1984；Nocedal和Wright，1999）：本发明方法让步骤1和步骤5运行直到满足最优条件（16）和（19）；接着，在步骤4，根据来自以前非约束内循环迭代的 是否满足以下两个条件中的一个，通过移除或添加一个或更多（网格）指标，更新有效集的总数：
[0093]
[0094]
[0095] 相反，在常规有效集方法中，一次仅更新有效集中的一个指标。换句话说，近似解一碰到上边界或下边界，常规有效集方法就更新有效集，但是本发明方法直到算法收敛到允许误差才更新有效集。常规有效集算法确保余量是单调递减的；在实际应用中，其示出非常缓慢的收敛速率，这是因为其往往退化到重新启动的最速下降法，尤其是在矩阵是高度病态时。这个新类型的有效集算法不保证给出单调递减数据失配。不过，经验表明，如果选择恰当允许误差，则这个算法在快速收敛速率方面完成的非常好并且产生满足允许误差的可行解。
[0096] 如同上述步骤1-6所例示的，本发明有效集类型算法的另一个优势是用于启动优化的开始的Q模型可以是任意的，这意味着开始的Q模型可以在上或下边界或甚至在边界之外设置。这个特征在某些情况下是非常有用的。
[0097] 本发明可以根据图3所示流程图实施。如同在常规的基于射线的层析成像算法中一样，在步骤60，速度模型20和源/接收器位置30被输入到射线追迹代码中，以输出射线路径信息，其用于在步骤70中构建核心矩阵（在方程式10和11中的A）。在步骤40，地震数据被分析并且源振幅谱由加权频率指数函数（方程式3）近似，在此期间对称指标n和特征频率f0被确定。接着在步骤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的功能。
[0098] 可以添加各种可选步骤或改进到图3所示的基本步骤。例如：
[0099] 1）所有地震轨迹可以在步骤40之前被预处理：第一批到达者被分窗口并被隔离，计算隔离子波的振幅谱，以及在振幅谱上实施多项式拟合。（原始地震数据往往是噪声非常多的，因此，不是利用高斯函数尝试拟合原始地震振幅谱，在实践中的第一步可以是原始地震振幅谱的多项式拟合，以将其转换为平滑曲线。在应用本发明方法中，这个相同中间步骤可能是有用的，这之后，用加权频率指数函数拟合平滑多项式曲线就更加容易。）[0100] 2）可以添加正则矩阵，以稳定步骤100的优化过程。
[0101] 3）可以从一个或更多速度异常建造开始的Q模型（步骤80）。
[0102] 4）当地震轨迹数量太小时，可以采用数据规则化技术。
[0103] 5）在步骤100期间，可以利用共轭梯度方法或LSQR方法（共轭梯度方法的变种；参见Paige和Saunders，1982）实施约束优化算法的内迭代。
[0104] 6）在步骤40期间，如果源谱不可获得，那么可以从接收的地震波形估算它，接收的地震波形的关联射线不穿越带有可能的Q异常的区域。
[0105] 示例
[0106] 在这部分，呈现Q层析成像的合成示例。本发明的Q层析成像算法的性能与利用带有或不带有箱约束的高斯函数的常规形心频率迁移Q层析成像算法的性能作比较。图4示出不对称的源振幅谱。源振幅谱的形心频率可以被数值地计算为55.7Hz，而峰值频率是2
33Hz。通过高斯函数源振幅谱拟合，高斯函数的峰值频率是55.7Hz而方差是1246Hz。通过加权频率指数函数拟合，在方程式（3）中的参数n和f0分别是1.45和22.7Hz。图5A是速度模型和通过运行射线追迹代码获得的射线路径，而图5B是真实的1/Q模型。在这个合成数据测试中，使用了50个源。对于每个源，存在90个接收器。因此，总共有4500个轨迹。
图5A-B到8A-B是其中绘出量的显示黑白复制品，根据绘图中的颜色标度，用颜色表示的深度和距离。
33Hz。通过高斯函数源振幅谱拟合，高斯函数的峰值频率是55.7Hz而方差是1246Hz。通过加权频率指数函数拟合，在方程式（3）中的参数n和f0分别是1.45和22.7Hz。图5A是速度模型和通过运行射线追迹代码获得的射线路径，而图5B是真实的1/Q模型。在这个合成数据测试中，使用了50个源。对于每个源，存在90个接收器。因此，总共有4500个轨迹。
图5A-B到8A-B是其中绘出量的显示黑白复制品，根据绘图中的颜色标度，用颜色表示的深度和距离。
[0107] 图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值以及存在很少人工产物。
[0108] 在本发明的某些实施例中，通过计算不穿越示出Q异常的地下区域的地震轨迹的平均振幅谱，执行对源振幅谱的估算。平滑可以应用于平均振幅谱，其可以利用多项式拟合技术完成。
[0109] 在本发明的某些实施例中，可以从地震数据轨迹的第一批到达者计算接收信号的振幅谱，这些地震数据轨迹可以通过窗口被隔离。
[0110] 在本发明的某些实施例中，可以从地下声学速度模型建造开始的1/Q模型。可以通过线性映射从地下声学速度模型中的一个或更多异常建造开始的1/Q模型。
[0111] 在本发明的某些实施例中，1/Q范围地图可以具有都大于0的下边界。
[0112] 前面的描述以本发明特定实施例为例子进行说明，不过对于本领域的技术人员显而易见的是，针对本文所述实施例的许多修改和变化是可能的。正如所附权利要求定义的一样，所有这样的修改和变化都应当在本发明的范围内。
[0113] 参考文献
[0114] Bierlaire,M.,Ph.L.Toint, 和 D.Tuyttens,“On iterative algorithms for linear least squares problems with bound constraints,”Linear Algebra Appl143,111-143(1991).
[0115] Delbos,D.,J.Ch.Gilbert,R.Glowinski, 和 D.Sinoquet,“Constrained optimization in seismic reflection tomography:a Gauss-Newton augmented Lagrangian approach,”Geophysics 164,670-684(2006).
[0116] Gao,F.,G.Fradelizio,A.Levander,G.Pratt, 和 C.Zelt,“Seismic velocity,Q,geological structure and lithology estimation at a ground water contamination site,”75th SEG，Expanded Abstracts,1561-1564,Soc.of Expl.Geophys(2005).
[0117] Hicks,G.J. 和 R.G.Pratt,“Reflection waveform inversion using local descent methods:Estimating attenuation and velocity over a gas-sand deposit，”Geophysics 66,598–612(2001).
[0118] Liao,Q. 和 G.A.McMechan,“Multifrequency viscoacoustic modeling and inversion,”Geophysics 61,No.5,1371-1378(1996).
[0119] P.,“Solving the minimal least squares problem subject tobounds on the variables,”BIT 24(1984).
[0120] Morigi,S.,L.Reichel,F.Sgallari, 和 F.Zama,“An iterative method for linear discrete ill-posed problems with box constraints,”Journal of Computational and Applied Mathematics 198,505-520(2007).
[0121] Nocedal,J. 和 S.J.Wright,Numerical Optimization,Springer,New York(1999).
[0122] Paige,C.C. 和 M.A.Saunders,“LSQR:An algorithm for sparse linear equations and sparse least squares,”ACM Trans.Math.Software 8,43-71(1982).[0123] Plessix,R.-E.,“Estimation of velocity and attenuation coefficient maps from crosswell seismic data,”Geophysics 71,S235–S240(2006).
[0124] Pratt,R.G.,K.Bauer， 和 M.Weber,“Crosshole waveform tomography velocity and attenuation images of arctic gas hydrates,”73rd SEG，Expanded Abstracts,Society of Exploration Geophysics,2255-2258(2003).
[0125] Rickett,J.E.,“Method for estimation of interval seismic quality factor,”U.S.Patent No.7,376,517(2006).
[0126] Rossi,G.,D.Gei,G. G.Madrussani, 和 J.M.Carcione,“Attenuation tomography:An application to gas-hydrate and free-gas detection,”Geophysics55,655-669(2007).
[0127] Spencer，T.W.,J.R.Sonnad, 和 T.M.Butler，“Seismic Q-stratigraphy or dissipation,”Geophysics 47,16-24(1982).
[0128] Tonn,R.,“The determination of the seismic quality factor Q from VSP data:A comparison of different computational method,”Geophysical Prospecting39,1–27(1991).
[0129] Quan,Y. 和 J.M.Harris,“Seismic attenuation tomography using the frequency shift method,”Geophysics 62,895-905(1997).
[0130] Watanabe,T.,K.T.Nihei,S.Nakagawa, 和 L.R.Myer，“Viscoacoustic waveform inversion of transmission data for velocity and attenuation,”J.Acoust.Soc.Am.115,3059–3067(2004).
[0131] Wright,S.J.,“Primal-Dual Interior Point Methods,”SIAM(1997).[0132] Wu,H., 和 J.M.Lees,“Attenuation structure of Coso geothermal data,California from wave pulse widths,”Bulletin of the Seismological Society of America 86,1574-1590(1996).
[0133] Zhang,C.,“Seismic Absorption Estimation and Compensation,”PhD thesis,University of British Columbia(2008).
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 光学相干层析成像装置 | 2020-05-12 | 327 |
| 一种光干涉层析成像仪 | 2020-05-12 | 559 |
| 光学相干层析成像系统 | 2020-05-12 | 719 |
| 医学热层析成像系统 | 2020-05-13 | 537 |
| 软场层析成像系统及软场层析成像方法 | 2020-05-13 | 772 |
| 光学层析成像装置 | 2020-05-12 | 523 |
| 光学层析成像装置 | 2020-05-12 | 327 |
| 中子相衬层析成像装置 | 2020-05-13 | 585 |
| 一种电容层析成像系统 | 2020-05-13 | 770 |
| 中子衍射层析成像装置 | 2020-05-13 | 585 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
QQ群二维码
意见反馈
意见反馈