多分量地震资料保幅角度域共成像点道集提取方法及系统与流程

1.本发明涉及多分量地震资料叠前深度偏移成像处理领域，特别是涉及一种 多分量地震资料保幅角度域共成像点道集提取方法及系统。


背景技术：

2.目前，全球油气勘探对象已逐步由常规向非常规、中深层向深层、中浅海 向深海、常规地带向极端地带转变。新的对象面临更深、更远、更极端环境等 挑战，是未来勘探技术攻关的重点和难点。地震偏移成像技术是地震勘探在油 气田开发应用中的关键技术之一。面对新的勘探对象，探测目标已由传统构造 油气藏转向隐蔽岩性油气藏，适用于传统构造油气藏的纵波偏移成像技术因利 用波场信息单一，在隐蔽岩性油气藏的成像应用中遇到瓶颈。多波偏移成像技 术利用了更丰富的地震波场信息，在小断层、缝洞、气云区及薄层等复杂构造 成像中显示出独特的优势和巨大的应用潜力。
3.共成像点道集提取是地震偏移成像技术的关键组成部分。共成像点道集所 包含的反映地下介质速度和岩性变化的信息是进行偏移速度分析和地震振幅 分析的有力工具。在偏移速度分析中，共成像点道集的剩余时差有效地反映了 偏移速度场的误差，当速度场准确时，共成像点道集同相轴是水平的，当速度 场不准确时，共成像点道集同相轴是弯曲的，而且弯曲的程度与速度场误差量 是有关系的，通过这种关系并借助共成像点道集拉平原则进行速度场的更新， 最终实现偏移速度建模。在地震振幅分析中，共成像点道集是叠前地震数据和 地震属性解释之间的桥梁。作为岩性分析的基础资料，共成像点道集的质量直 接决定了地震振幅分析技术能否准确的分析和描述岩性油气藏，也决定了弹性 参数及储层参数反演的成败。面对当前复杂的油气勘探环境，传统的偏移距域 共成像点道集因多路径问题会产生严重假象，在偏移速度分析和地震振幅分析 方面有较大的局限性，为了解决这一问题，学者们提出了角度域共成像点道集。
4.角度域共成像点道集因不受多路径问题影响，是一种具有良好应用潜力的 共成像点道集。基于逆时偏移的角度域共成像点道集因其采用精确的地震波动 理论，能够适用于任意复杂构造，且算法精度高、稳健性强，成为了当前油气 勘探界应用的主流。然而，由于受速度不准、采集孔径有限、照明不均以及偏 移算子不精确等因素影响，基于当前已有方法所提取的角度域共成像点道集分 辨率低、振幅不均衡、低波数假象严重，这对后续的数据处理和解释造成很大 影响，使其很难直接用于实际生产，也使得角度域共成像点道集的优势难以发 挥。为此，必须针对多分量地震资料建立一套新的、可提取高质量的多波角度 域共成像点道集的成像方法及系统。


技术实现要素：

5.本发明的目的是提供一种多分量地震资料保幅角度域共成像点道集提取 方法及系统，可以获得高精度、高分辨率、高信噪比、振幅保真的多波角度域 共成像点道集，为地震叠前岩性反演及流体识别与预测提供高质量基础数据， 指导勘探部署。
6.为实现上述目的，本发明提供了如下方案：
7.一种多分量地震资料保幅角度域共成像点道集提取方法，包括：
8.依据勘探目标及工区的地质地球物理条件，确定进行角度域共成像点道集 提取所需的观测的多炮多分量观测地震记录、纵波偏移速度模型、横波偏移速 度、偏移密度模型及观测系统参数；
9.针对第k炮，根据所述纵波偏移速度模型、所述横波偏移速度模型、所述 偏移密度模型及所述观测系统参数，基于数值方法求解地震波动方程，构建正 传多分量地震波场；
10.针对第k炮，根据所述正传多分量地震波场和所述多分量观测地震记录基 于数值方法求解地震波场伴随方程，计算第i次迭代的偏移距域共成像点道集 梯度剖面；
11.基于所述第i次迭代的多波局部偏移距域共成像点道集梯度剖面，构建第 i次迭代的偏移距域共成像点道集下降方向剖面；
12.利用所述下降方向剖面，基于数值方法求解以多波偏移距域共成像点道集 为散射源的波场正传方程，得到预测纵波地震波场和预测横波地震波场；
13.基于所述预测纵波地震波场和所述预测横波地震波场，利用波场合成算子 计算预测多分量地震记录增量；
14.估算最优化步长；
15.根据所述最优化步长、所述下降方向剖面和所述多分量地震记录增量，更 新偏移距域共成像点道集和预测多分量地震记录；
16.判定是否满足收敛标准；
17.若是，则输出更新后的偏移距域共成像点道集，所述更新后的偏移距域共 成像点道集为最终偏移距域共成像点道集；
18.根据所述最终偏移距域共成像点道集，确定最优角度域共成像点道集；
19.若否，则针对第k炮，根据所述多分量背景地震波场和所述多分量观测地 震记录基于数值方法求解地震波场伴随方程，计算第i+1次迭代的偏移距域共 成像点道集梯度剖面，直至获得最优角度域共成像点道集。
20.可选地，所述针对第k炮，根据所述纵波偏移速度模型、所述横波偏移速 度模型、所述偏移密度模型及所述观测系统参数，基于数值方法求解地震波动 方程，构建正传多分量地震波场，具体包括：
21.针对第k炮，基于所述观测系统参数，得到所述第k炮的炮点坐标，在所 述第k炮对应炮点位置设置震源子波，利用所述纵波偏移速度模型、所述横波 偏移速度模型、所述偏移密度模型及所述观测系统参数，基于数值方法求解地 震波动方程，得到第k炮的正传多分量地震波场。
22.可选地，所述针对第k炮，根据所述正传多分量地震波场和所述多分量观 测地震记录基于数值方法求解地震波场伴随方程，计算第i次迭代的偏移距域 共成像点道集梯度，具体包括：
23.设置当前的迭代次数i，针对第k炮，利用当前炮的观测多分量地震记录 及第i
‑
1次迭代更新所得的预测多分量地震记录，得到多分量地震记录残差；
24.以所述多分量地震记录残差为边值条件，基于数值方法逆时求解地震纯纵 波波动方程，得到第k炮的反传多分量地震波场；
25.在相同的时刻，对所述第k炮的正传多分量地震波场和所述第k炮的反传 多分量地震波场进行地震波模式分解，获得相同时刻的标量纵波正传地震波 场、标量横波正传地震波场、标量纵波反传地震波场和标量横波反传地震波场；
26.对所述标量纵波正传地震波场、所述标量横波正传地震波场、所述标量纵 波反传地震波场和所述标量横波反传地震波场应用多波局部偏移距域共成像 点道集梯度计算方程，得到单炮多波局部偏移距域共成像点道集梯度；
27.根据所有所述单炮多波局部偏移距域共成像点道集梯度，依据观测系统位 置信息进行叠加，得到第i次迭代的多波局部偏移距域共成像点道集梯度剖面。
28.可选地，所述基于所述第i次迭代的多波局部偏移距域共成像点道集梯度 剖面，构建第i次迭代的偏移距域共成像点道集下降方向剖面，具体包括：
29.基于第i次迭代的多波局部偏移距域共成像点道集梯度剖面，利用反演方 法，得到第i次迭代的偏移距域共成像点道集下降方向剖面。
30.可选地，所述利用所述下降方向剖面，基于数值方法求解以多波偏移距域 共成像点道集为散射源的波场正传方程，得到预测纵波地震波场和预测横波地 震波场，具体包括：
31.利用第k炮的观测系统位置信息，从第i次迭代的偏移距域共成像点道集 下降方向剖面中读取相同位置处的下降方向；
32.采用地震波模式分解方程对所述正传多分量地震波场进行波模式分解，得 到第k炮的标量纵波正传地震波场和第k炮的标量横波正传地震波场；
33.基于第k炮的纵波下降方向和第k炮的标量纵波正传地震波场，确定第k 炮的预测纵波地震波场；
34.基于第k炮的横波下降方向和第k炮的标量横波正传地震波场，确定第k 炮的预测横波地震波场。
35.可选地，所述根据所述最优化步长、所述下降方向剖面和所述多分量地震 记录增量，更新偏移距域共成像点道集和预测多分量地震记录，具体包括：
36.根据所述最优化步长和所述下降方向剖面，更新第i次迭代的偏移距域共 成像点道集；
37.根据所述最优化步长和所述多分量地震记录增量，更新第i次迭代的预测 多分量地震记录。
38.可选地，所述收敛标准具体为：
[0039][0040]
其中，relerr为迭代停止的阈值标准，relerr选取1.0
e
‑4，misfit
i
为第i次迭代 的目标函数值，misfit
i
‑1为第i
‑
1次迭代的目标函数值，d
i
为 第k炮第i次迭代的预测地震记录，d
w
为第k炮的观测纯波地震记录。
[0041]
一种多分量地震资料保幅角度域共成像点道集提取系统，包括：
[0042]
数据获取模块，用于依据勘探目标及工区的地质地球物理条件，确定进行 角度域共成像点道集提取所需的观测的多炮多分量观测地震记录、纵波偏移速 度模型、横波偏移
速度、偏移密度模型及观测系统参数；
[0043]
正传多分量地震波场构建模块，用于针对第k炮，根据所述纵波偏移速度 模型、所述横波偏移速度模型、所述偏移密度模型及所述观测系统参数，基于 数值方法求解地震波动方程，构建正传多分量地震波场；
[0044]
梯度剖面第一确定模块，用于针对第k炮，根据所述正传多分量地震波场 和所述多分量观测地震记录基于数值方法求解地震波场伴随方程，计算第i次 迭代的偏移距域共成像点道集梯度剖面；
[0045]
下降方向剖面确定模块，用于基于所述第i次迭代的多波局部偏移距域共 成像点道集梯度剖面，构建第i次迭代的偏移距域共成像点道集下降方向剖面；
[0046]
预测纵/横波地震波场确定模块，用于利用所述下降方向剖面，基于数值 方法求解以多波偏移距域共成像点道集为散射源的波场正传方程，得到预测纵 波地震波场和预测横波地震波场；
[0047]
预测多分量地震记录增量确定模块，用于基于所述预测纵波地震波场和所 述预测横波地震波场，利用波场合成算子计算预测多分量地震记录增量；
[0048]
最优化步长估算模块，用于估算最优化步长；
[0049]
更新模块，用于根据所述最优化步长、所述下降方向剖面和所述多分量地 震记录增量，更新偏移距域共成像点道集和预测多分量地震记录；
[0050]
判断模块，用于判定是否满足收敛标准；
[0051]
最终偏移距域共成像点道集输出模块，用于满足收敛标准时，则输出更新 后的偏移距域共成像点道集，所述更新后的偏移距域共成像点道集为最终偏移 距域共成像点道集；
[0052]
最优角度域共成像点道集确定模块，用于根据所述最终偏移距域共成像点 道集，确定最优角度域共成像点道集；
[0053]
梯度剖面第二确定模块，用于在不满足收敛标准时，针对第k炮，根据所 述多分量背景地震波场和所述多分量观测地震记录基于数值方法求解地震波 场伴随方程，计算第i+1次迭代的偏移距域共成像点道集梯度剖面。
[0054]
根据本发明提供的具体实施例，本发明公开了以下技术效果：
[0055]
1)本发明为一种多分量地震资料保幅角度域共成像点道集提取方法，与 常规角度域共成像点道集提取方法相比，本发明可以获得高精度、高分辨率、 高信噪比、振幅保真的多波角度域共成像点道集；2)本发明以多波局部偏移 距域共成像点道集为反演目标，通过弹性波逆时偏移和弹性波逆时反偏移，可 以获得高质量的多波局部偏移距域共成像点道集，再基于局部偏移距域共成像 点道集与角度域共成像点道集的转换关系，获得多波角度域共成像点道集，该 多波角度域共成像点道集，可以直接用于后续的地震资料叠前反演，极大提高 了储层预测及流体识别的精度。
附图说明
[0056]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是 本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的 前提下，还可以根据这些附
图获得其他的附图。
[0057]
图1为本发明多分量地震资料保幅角度域共成像点道集提取方法流程图；
[0058]
图2是本发明提供的二维水平层状介质模型模型；
[0059]
图3是图2所示二维水平层状介质模型的多波局部偏移距域共成像点道 集：其中，图3(a)是利用传统方法所得的纵波局部偏移距域共成像点道集， 图3(b)是利用传统方法所得的横波局部偏移距域共成像点道集，图3(c) 是利用本发明所得的纵波局部偏移距域共成像点道集，图3(d)是利用本发 明所得的横波局部偏移距域共成像点道集；
[0060]
图4是图2所示二维水平层状介质模型的多波角度域共成像点道集：其中， 图4(a)是利用传统方法所得的纵波角度域共成像点道集，图4(b)是利用 传统方法所得的横波角度域共成像点道集，图4(c)是利用本发明所得的纵 波角度域共成像点道集，图4(d)是利用本发明所得的横波角度域共成像点 道集；
[0061]
图5是本发明提供的二维简单介质模型；
[0062]
图6是图5所示的二维简单介质模型的地震记录剖面；其中，图6(a) 第20炮的观测记录x分量，图6(b)第20炮的观测记录z分量，图6(c) 利用传统方法所得的第20炮的预测记录x分量，图6(d)利用传统方法所得 的第20炮的预测记录z分量，图6(e)利用本发明所得的第20炮的预测记 录x分量，图6(f)利用本发明所得的第20炮的预测记录z分量；
[0063]
图7是图5所示的二维简单介质模型的多波局部偏移距域共成像点道集； 其中，图7(a)是利用传统方法所得的纵波局部偏移距域共成像点道集，图7 (b)是利用传统方法所得的横波局部偏移距域共成像点道集，图7(c)是利 用本发明所得的纵波局部偏移距域共成像点道集，图7(d)是利用本发明所 得的横波局部偏移距域共成像点道集；
[0064]
图8是图5所示的二维简单介质模型的多波角度域共成像点道集；其中， 图8(a)是利用传统方法所得的纵波角度域共成像点道集，图8(b)是利用 传统方法所得的横波角度域共成像点道集，图8(c)是利用本发明所得的纵 波角度域共成像点道集，图8(d)是利用本发明所得的横波角度域共成像点 道集；
[0065]
图9是本发明提供的marmousi
‑
2偏移模型；其中，(a)纵波速度模型， (b)横波速度模型，(c)密度模型；
[0066]
图10是图9所示marmousi
‑
2模型的多炮叠加偏移剖面：其中，图10(a) 第55炮的观测记录x分量，图10(b)第55炮的观测记录z分量，图10(c) 利用传统方法所得的第55炮的预测记录x分量，图10(d)利用传统方法所 得的第55炮的预测记录z分量，图10(e)利用本发明所得的第55炮的预测 记录x分量，图10(f)利用本发明所得的第55炮的预测记录z分量；
[0067]
图11是图9所示marmousi
‑
2模型的纵波局部偏移距域共成像点道集；其 中，图11(a)是利用传统方法所得的纵波局部偏移距域共成像点道集，图11 (b)是利用本发明所得的纵波局部偏移距域共成像点道集；
[0068]
图12是图9所示marmousi
‑
2模型的横波局部偏移距域共成像点道集；其 中，图12(a)是利用传统方法所得的横波局部偏移距域共成像点道集，图12 (b)是利用本发明所得的横波局部偏移距域共成像点道集；
[0069]
图13是图9所示marmousi
‑
2模型的纵波角度域共成像点道集；其中，图 13(a)是利用传统方法所得的纵波角度域共成像点道集，图13(b)是利用 本发明所得的纵波角度域共成像点道集；
[0070]
图14是图9所示marmousi
‑
2模型的横波角度域共成像点道集；其中，图 14(a)是利用传统方法所得的横波角度域共成像点道集，图14(b)是利用 本发明所得的横波角度域共成像点道集；
[0071]
图15是图9所示marmousi
‑
2模型的角度域共成像点道集叠加偏移剖面： 其中，图15(a)是利用传统方法所得的纵波角度域共成像点道集叠加偏移剖 面，图15(b)利用传统方法所得的横波角度域共成像点道集叠加偏移剖面， 图15(c)是利用本发明所得的纵波角度域共成像点道集叠加偏移剖面，图15 (d)是利用本发明所得的横波角度域共成像点道集叠加偏移剖面。
具体实施方式
[0072]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清 楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是 全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造 性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0073]
本发明的目的是提供一种多分量地震资料保幅角度域共成像点道集提取 方法及系统，可以获得高精度、高分辨率、高信噪比、振幅保真的多波角度域 共成像点道集，为地震叠前岩性反演及流体识别与预测提供高质量基础数据， 指导勘探部署。
[0074]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和 具体实施方式对本发明作进一步详细的说明。
[0075]
图1为本发明多分量地震资料保幅角度域共成像点道集提取方法流程图。 如图1所示，一种多分量地震资料保幅角度域共成像点道集提取方法包括：
[0076]
步骤101：依据勘探目标及工区的地质地球物理条件，确定进行角度域共 成像点道集提取所需的观测的多炮多分量观测地震记录、纵波偏移速度模型、 横波偏移速度、偏移密度模型及观测系统参数。
[0077]
具体的，依据勘探目标及工区的地质地球物理条件，确定进行角度域共成 像点道集提取所需的观测的多炮多分量观测地震记录d(x
r
,t；x
s
)＝(d
x
,d
y
,d
z
)、纵 波偏移速度模型v
p
(x)、横波偏移速度v
s
(x)、偏移密度模型ρ(x)；确定观测系 统参数，建立目标靶区在笛卡尔坐标系下的观测系统；其中，d
x
、d
y
和d
z
分 别表示观测的多分量观测地震记录d在笛卡尔坐标系下x、y和z方向的分量， x＝(x,y,z)表示地下网格位置坐标矢量，x
s
＝(x
s
,y
s
,z
s
)表示震源空间位置矢量， x
r
＝(x
r
,y
r
,z
r
)表示检波点空间位置矢量，t表示波传播时间。
[0078]
步骤102：针对第k炮，根据所述纵波偏移速度模型、所述横波偏移速度 模型、所述偏移密度模型及所述观测系统参数，基于数值方法求解地震波动方 程，构建正传多分量地震波场，具体包括：
[0079]
针对第k炮，基于所述观测系统参数，得到所述第k炮的炮点坐标，在所 述第k炮对应炮点位置设置震源子波，利用所述纵波偏移速度模型、所述横波 偏移速度模型、所述偏移密度模型及所述观测系统参数，基于数值方法求解地 震波动方程，得到第k炮的正传多分量地震波场。
[0080]
具体的，针对第k炮，所述第k炮为任意一炮，基于确定的观测系统参数 获得该炮的炮点坐标，在该炮对应炮点位置设置震源子波，利用上述确定的纵 波偏移速度模型、横
波偏移速度模型、偏移密度模型及观测系统参数，基于数 值方法求解地震波动方程，实现对该炮点波场的正演模拟，获得该炮的每一个 时刻的正传多分量地震波场进而获得该炮的每一个时刻 的正传多分量地震波场；其中，和分别表示多分量地震波场u
f
在笛 卡尔坐标系下x、y和z方向的分量。
[0081]
其中所述的基于数值方法求解地震波动方程，实现对该炮点波场的正演模 拟，获得该炮的每一个时刻的正传多分量地震波场具体 为利用偏微分方程数值方法求解：
[0082][0083]
所述方程(1)中，f(x,t；x
s
)＝(f
x
,f
y
,f
z
)表示震源矢量；x、y和z分别表示笛 卡尔坐标系下的x方向、y方向和z方向；参数λ(x)和μ(x)具体为：
[0084][0085][0086]
步骤103：针对第k炮，根据所述正传多分量地震波场和所述多分量观测 地震记录基于数值方法求解地震波场伴随方程，计算第i次迭代的偏移距域共 成像点道集梯度剖面，具体包括：
[0087]
步骤1031：设置当前的迭代次数i，针对第k炮，利用当前炮的观测多分 量地震记录及第i
‑
1次迭代更新所得的预测多分量地震记录，得到多分量地震 记录残差。
[0088]
步骤1032：以所述多分量地震记录残差为边值条件，基于数值方法逆时 求解地震纯纵波波动方程，得到第k炮的反传多分量地震波场。
[0089]
步骤1033：在相同的时刻，对所述第k炮的正传多分量地震波场和所述 第k炮的反传多分量地震波场进行地震波模式分解，获得相同时刻的标量纵波 正传地震波场、标量横波正传地震波场、标量纵波反传地震波场和标量横波反 传地震波场。
[0090]
步骤1034：对所述标量纵波正传地震波场、所述标量横波正传地震波场、 所述标量纵波反传地震波场和所述标量横波反传地震波场应用多波局部偏移 距域共成像点道集梯度计算方程，得到单炮多波局部偏移距域共成像点道集梯 度。
[0091]
步骤1035：根据所有所述单炮多波局部偏移距域共成像点道集梯度，依 据观测系统位置信息进行叠加，得到第i次迭代的多波局部偏移距域共成像点 道集梯度剖面。
[0092]
具体的，设置当前的迭代次数i，针对每一炮，利用该炮的观测多分量地 震记录d及第i
‑
1次迭代更新所得的预测多分量地震记录d
i
(x
r
,t；x
s
)＝(d
x
,d
y
,d
z
)， 计算多分量地
震记录残差δd
i
(x
r
,t；x
s
)＝(δd
x
,δd
y
,δd
z
)，其计算方法为 δd
i
(x
r
,t；x
s
)＝d
i
‑
d；以该炮的多分量地震记录残差δd
i
为边值条件，基于数值方 法逆时求解地震波动方程，获得该炮的每一个时刻的反传多分量地震波场 读取步骤102中计算的该炮的正传多分量地震波场，在 相同的时刻，对该炮的正传多分量地震波场和反传多分量地震波场进行地震波 模式分解，获得相同时刻的标量纵波正传地震波场p
f
(x,t；x
s
)、标量横波正传地 震波场s
f
(x,t；x
s
)、标量纵波反传地震波场p
b
(x,t；x
s
)和标量横波反传地震波场 s
b
(x,t；x
s
)；对该炮的所述标量纵波正传地震波场、标量横波正传地震波场、标 量纵波反传地震波场和标量横波反传地震波场应用多波局部偏移距域共成像 点道集梯度计算方程，获得该炮的单炮多波局部偏移距域共成像点道集梯度， 进而获得每一炮对应的多波局部偏移距域共成像点道集梯度；所有单炮多波局 部偏移距域共成像点道集依据观测系统位置信息进行叠加构成了本次迭代的 多波局部偏移距域共成像点道集梯度具体包括纵波局部偏移 距域共成像点道集梯度和横波局部偏移距域共成像点道集梯度其中，d
x
、d
y
和d
z
分别表示预测的多分量地震记录d在笛卡尔坐标系下x、y 和z方向的分量，δd
x
、δd
y
和δd
z
分别表示多分量地震记录残差δd在笛卡尔 坐标系下x、y和z方向的分量，和分别表示反传多分量地震波场u
b
在笛卡尔坐标系下x、y和z方向的分量，h＝(h
x
,h
y
,h
z
)表示局部偏移距矢量， h
x
、h
y
和h
z
分别表示局部偏移距矢量h在笛卡尔坐标系下的x方向、y方向和 z方向分量；
[0093]
所述的以该炮的多分量地震记录残差δd
i
为边值条件，基于数值方法逆时 求解地震波动方程，获得该炮的每一个时刻的反传多分量地震波场 具体为利用偏微分方程数值方法求解：
[0094][0095]
所述的在相同的时刻，对该炮的正传多分量地震波场和反传多分量地震波 场进行地震波模式分解，获得相同时刻的标量纵波正传地震波场p
f
(x,t；x
s
)、标 量横波正传地震波场s
f
(x,t；x
s
)、标量纵波反传地震波场p
b
(x,t；x
s
)和标量横波反 传地震波场s
b
(x,t；x
s
)，具体为：
[0096]
[0097]
和
[0098][0099]
所述方程(5)和(6)中，f
‑1{
·
}表示傅里叶逆变换算子；a
p
(k,t；x
s
)和a
s
(k,t；x
s
) 分别表示波数域的纵波偏振矢量和横波偏振矢量；k＝(k
x
,k
y
,k
z
)表示波矢量； 和分别表示波数域的正传多分量地震波场和反传多分量地 震波场；“·”表示矢量内积运算；所述波数域的纵波偏振矢量和横波偏振矢量， 具体为求解：
[0100][0101]
所述方程(7)中，角标b表示波模式，包含纵波和横波；g是克里斯托 弗矩阵，其元素g
jk
＝c
jklm
n
k
n
m
，c
jklm
是介质刚度矩阵，n为归一化波数，下角 标的字母j、k、l和m均可以取x、y和z，表示笛卡尔坐标系下x、y和z方 向；
[0102]
所述的对该炮的所述标量纵波正传地震波场、标量横波正传地震波场、标 量纵波反传地震波场和标量横波反传地震波场应用多波局部偏移距域共成像 点道集梯度计算方程，获得该炮的单炮多波局部偏移距域共成像点道集梯度， 具体为：
[0103][0104]
所述方程(8)即为多波局部偏移距域共成像点道集梯度计算方程。
[0105]
步骤104：基于所述第i次迭代的多波局部偏移距域共成像点道集梯度剖 面，构建第i次迭代的偏移距域共成像点道集下降方向剖面，具体包括：
[0106]
基于第i次迭代的多波局部偏移距域共成像点道集梯度剖面，利用反演方 法，得到第i次迭代的偏移距域共成像点道集下降方向剖面其中，表示第i次迭代的纵波下降方向，表示第i次迭代的横波 下降方向。
[0107]
步骤105：利用所述下降方向剖面，基于数值方法求解以多波偏移距域共 成像点道集为散射源的波场正传方程，得到预测纵波地震波场和预测横波地震 波场，具体包括：
[0108]
步骤1051：利用第k炮的观测系统位置信息，从第i次迭代的偏移距域共 成像点道集下降方向剖面中读取相同位置处的下降方向。
[0109]
利用第k炮的观测系统位置信息，从步骤104中所得的第i次迭代的下降 方向dg
i
(x,h)中读取相同位置处的该炮的下降方向dg(x,h；x
s
)，具体包括该炮的 纵波下降方向dg
p
(x,h；x
s
)和横波下降方向dg
s
(x,h；x
s
)。
[0110]
步骤1052：采用地震波模式分解方程对所述正传多分量地震波场进行波 模式分解，得到第k炮的标量纵波正传地震波场和第k炮的标量横波正传地震 波场。即：
[0111]
读取步骤102中计算的该炮的正传多分量地震波场，利用方程(5)所示 的地震波模式分解方程对该炮的正传多分量地震波场进行波模式分解，获得该 炮的标量纵波正传地震波场p
f
(x,t；x
s
)和标量横波正传地震波场s
f
(x,t；x
s
)。
[0112]
步骤1053：基于第k炮的纵波下降方向和第k炮的标量纵波正传地震波 场，确定第k炮的预测纵波地震波场。即：
[0113]
基于该炮的纵波下降方向dg
p
(x,h；x
s
)和该炮的标量纵波正传地震波场 p
f
(x,t；x
s
)，利用上述确定的纵波偏移速度模型及偏移参数，基于数值方法，应 用纵波散射条件，求解纵波地震散射方程，实现对该炮点预测纵波地震波场的 数值模拟，获得该炮的每一个时刻的预测纵波地震波场p
c
(x,t；x
s
)。其中纵波散 射条件具体为：
[0114][0115]
所述方程(9)中，f
p
(x,t；x
s
)表示纵波散射源；基于所述方程(9)表示的 纵波散射源，所述的纵波地震散射方程具体为：
[0116][0117]
步骤1054：基于第k炮的横波下降方向和第k炮的标量横波正传地震波 场，确定第k炮的预测横波地震波场。即：
[0118]
基于该炮的横波下降方向dg
s
(x,h；x
s
)和该炮的标量横波正传地震波场 s
f
(x,t；x
s
)，利用上述确定的横波偏移速度模型及偏移参数，基于数值方法，应 用横波散射条件，求解横波地震散射方程，实现对该炮点预测横波地震波场的 数值模拟，获得该炮的每一个时刻的预测横波地震波场s
c
(x,t；x
s
)。其中横波散 射条件具体为：
[0119][0120]
所述方程(11)中，f
s
(x,t；x
s
)表示横波散射源；基于所述方程(11)表示 的横波散射源，所述的横波地震散射方程具体为：
[0121][0122]
步骤106：基于所述预测纵波地震波场和所述预测横波地震波场，利用波 场合成算子计算预测多分量地震记录增量。即：
[0123]
在相同的时刻，利用多分量地震波场合成方程对所得的预测纵波地震波场 p
c
(x,t；x
s
)和预测横波地震波场s
c
(x,t；x
s
)进行多分量地震波场合成，获得该炮的 每一个时刻的预测多分量地震波场基于该炮的观测系统 信息，对该炮的预测多分量地震波场进行采样，获得第该炮的在第i次迭代的 多分量地震记录增量δd
i
(x
r
,t；x
s
)＝(δd
x
,δd
y
,δd
z
)；进而获得每一炮在第i次迭代 的多分量地震记录增量；其中，δd
x
、δd
y
和δd
z
分别表示多分量地震记录增量 δd在笛卡尔坐标系下x、y和z方向的分量。
[0124]
具体为：
[0125][0126]
所述方程(13)中，和分别表示波数域的预测预测纵波 地震波场和波数域的预测横波地震波场；
[0127]
所述的基于该炮的观测系统信息，对该炮的预测多分量地震波场进行采 样，获得第该炮的在第i次迭代的多分量地震记录增量 δd
i
(x
r
,t；x
s
)＝(δd
x
,δd
y
,δd
z
)，具体为：
[0128]
δd
i
(x
r
,t；x
s
)＝u
c
(x
r
,t；x
s
)；
ꢀꢀꢀꢀꢀ
(14)
[0129]
步骤107：估算最优化步长，所述最优化步长用α
i
表示。
[0130]
步骤108：根据所述最优化步长、所述下降方向剖面和所述多分量地震记 录增量，更新偏移距域共成像点道集和预测多分量地震记录，具体包括：
[0131]
步骤1081：根据所述最优化步长和所述下降方向剖面，更新第i次迭代的 偏移距域共成像点道集。即利用步骤107得到的最优化步长α
i
及步骤104得到 的下降方向剖面dg
i
，更新第i次迭代的多波局部偏移距域共成像点道集 r
i
(x,h)＝r
i
‑1(x,h)+α
i
dg
i
(x,h)；其中，多波局部偏移距域共成像点道集 包含纵波局部偏移距域共成像点道集和横波局部偏 移距域共成像点道集
[0132]
步骤1082：根据所述最优化步长和所述多分量地震记录增量，更新第i 次迭代的预测多分量地震记录。即利用步骤107得到的最优化步长α
i
及步骤 105得到的多分量地震记录增量δd
i
，更新第i次迭代的预测多分量地震记录 d
i
＝d
i
‑1+α
i
δd
i
，其中，d0＝0。
[0133]
步骤109：判定是否满足收敛标准。具体的，基于l
‑
2范数，计算第i次 迭代的目标函数值misfit
i
，判断当前迭代是否满足收敛标准。
[0134]
所述收敛标准具体为：
[0135][0136]
其中，relerr为迭代停止的阈值标准，relerr选取1.0
e
‑4，misfit
i
为第i次迭代的 目标函数值，misfit
i
‑1为第i
‑
1次迭代的目标函数值，d
i
为第k 炮第i次迭代的预测地震记录，d
w
为第k炮的观测纯波地震记录。
[0137]
步骤110：若是，则输出更新后的偏移距域共成像点道集，所述更新后的 偏移距域共成像点道集为最终偏移距域共成像点道集r(x,h)。
[0138]
步骤111：根据所述最终偏移距域共成像点道集，确定最优角度域共成像 点道集。
[0139]
利用步骤109所得的最优多波局部偏移距域共成像道集r(x,h)，基于多波 局部偏移距域共成像点道集与多波角度域共成像点道集的转换关系，获得最优 多波角度域共成像点道集r(x,θ)；其中，表示θ反射角。
[0140]
所述的多波局部偏移距域共成像点道集与多波角度域共成像点道集的转 换关系，具体是利用局部偏移距矢量h和反射角θ之间的映射关系，对多波局 部偏移距域共成像点道集进行数据重排，获得多波角度域共成像点道集，具体 为：对纵波局部偏移距域共成像点道集r(x,h)的每一个水平空间位置进行共中 心点道集抽取，获得r(x
′
,h)，其中x
′
＝(m
x
,m
y
,z)，m
x
和m
y
分别表示笛卡尔坐 标系下x方向和y方向的共中心点坐标。
[0141]
其次，对所抽取的共成像点道集r(x
′
,h)进行六维傅里叶变换，获得相应的 多波波数域纵波共成像点道集其中表示空间波数， 表示笛卡尔坐标系下x方向共中心点波数，表示笛卡尔坐标系下y方向 共中心点波数，k
z
表示笛
卡尔坐标系z方向波数；表示局部偏 移距波数，表示笛卡尔坐标系下x方向局部偏移距波数，表示笛卡尔坐 标系下y方向局部偏移距波数，表示笛卡尔坐标系下z方向局部偏移距波 数。
[0142]
然后，应用地震波传播反射角θ与波数k
′
x
和k
h
之间的关系，对多波波数域 的共成像点道集进行逐一投影，最终对投影后的共成像点道集进行六 维傅里叶逆变换，获得多波角度域共成像点道集梯度r(x,θ)，其中，所述的地 震波传播反射角θ与波数k
′
x
和k
h
之间的关系具体为：
[0143][0144]
步骤112：若否，则针对第k炮，根据所述多分量背景地震波场和所述多 分量观测地震记录基于数值方法求解地震波场伴随方程，计算第i+1次迭代的 偏移距域共成像点道集梯度剖面，直至获得最优的多波局部偏移距域共成像点 道集。利用最优的多波局部偏移距域共成像点道集为地震叠前岩性反演及流体 识别与预测提供高质量基础数据，指导勘探部署。
[0145]
对应于本发明的一种多分量地震资料保幅角度域共成像点道集提取方法， 本发明还提供一种多分量地震资料保幅角度域共成像点道集提取系统，该系统 包括：
[0146]
数据获取模块，用于依据勘探目标及工区的地质地球物理条件，确定进行 角度域共成像点道集提取所需的观测的多炮多分量观测地震记录、纵波偏移速 度模型、横波偏移速度、偏移密度模型及观测系统参数；
[0147]
正传多分量地震波场构建模块，用于针对第k炮，根据所述纵波偏移速度 模型、所述横波偏移速度模型、所述偏移密度模型及所述观测系统参数，基于 数值方法求解地震波动方程，构建正传多分量地震波场；
[0148]
梯度剖面第一确定模块，用于针对第k炮，根据所述正传多分量地震波场 和所述多分量观测地震记录基于数值方法求解地震波场伴随方程，计算第i次 迭代的偏移距域共成像点道集梯度剖面；
[0149]
下降方向剖面确定模块，用于基于所述第i次迭代的多波局部偏移距域共 成像点道集梯度剖面，构建第i次迭代的偏移距域共成像点道集下降方向剖面；
[0150]
预测纵/横波地震波场确定模块，用于利用所述下降方向剖面，基于数值 方法求解以多波偏移距域共成像点道集为散射源的波场正传方程，得到预测纵 波地震波场和预测横波地震波场；
[0151]
预测多分量地震记录增量确定模块，用于基于所述预测纵波地震波场和所 述预测横波地震波场，利用波场合成算子计算预测多分量地震记录增量；
[0152]
最优化步长估算模块，用于估算最优化步长；
[0153]
更新模块，用于根据所述最优化步长、所述下降方向剖面和所述多分量地 震记录增量，更新偏移距域共成像点道集和预测多分量地震记录；
[0154]
判断模块，用于判定是否满足收敛标准；
[0155]
最终偏移距域共成像点道集输出模块，用于满足收敛标准时，则输出更新 后的偏移距域共成像点道集，所述更新后的偏移距域共成像点道集为最终偏移 距域共成像点道集；
[0156]
最优角度域共成像点道集确定模块，用于根据所述最终偏移距域共成像点 道集，确定最优角度域共成像点道集；
[0157]
梯度剖面第二确定模块，用于在不满足收敛标准时，针对第k炮，根据所 述多分量背景地震波场和所述多分量观测地震记录基于数值方法求解地震波 场伴随方程，计算第i+1次迭代的偏移距域共成像点道集梯度剖面。
[0158]
实例1：
[0159]
图2是二维水平层状介质模型模型。在此模型上设置21个爆炸震源，震 源子波设定为雷克子波，主频为20赫兹，起始震源点的位于(2500m，100m) 处，炮间隔为50m。采用中间放炮两边接收观测系统，单边最大偏移距2500m， 最小偏移距为0m，道间距为10m，记录时间长度2.5s，时间步长1ms。图3 是图2所示二维水平层状介质模型的多波局部偏移距域共成像点道集：其中， 图3(a)是利用传统方法所得的纵波局部偏移距域共成像点道集，图3(b) 是利用传统方法所得的横波局部偏移距域共成像点道集，图3(c)是利用本 发明所得的纵波局部偏移距域共成像点道集，图3(d)是利用本发明所得的 横波局部偏移距域共成像点道集。由图3可知，传统方法所得多波局部偏移距 域共成像点道集存在明显的子波旁瓣，而且远离零偏移距处有较强的能量残 留。剖面存在比较明显的噪声，剖面的分辨率较低，振幅不均衡。相比之下， 本发明所得的多波局部偏移距域共成像点道集聚焦性更好，子波旁瓣得到压 制，假象更少，也证明了本发明的可行性及有效性。图4是图2所示二维水平 层状介质模型的多波角度域共成像点道集：其中，图4(a)是利用传统方法 所得的纵波角度域共成像点道集，图4(b)是利用传统方法所得的横波角度 域共成像点道集，图4(c)是利用本发明所得的纵波角度域共成像点道集， 图4(d)是利用本发明所得的横波角度域共成像点道集。由图4可以看出， 利用本发明所得的多波角度域共成像点道集具有更高的分辨率、信噪比及振幅 均衡性，这间接证明了本发明的有效性。
[0160]
实例2：
[0161]
图5是二维简单介质模型。在此模型上设置39个爆炸震源，震源子波设 定为雷克子波，主频为20赫兹，起始震源点的位于(100m，100m)处，炮 间隔为100m。采用中间放炮两边接收观测系统，单边最大偏移距1500m，最 小偏移距为0m，道间距为10m，记录时间长度2.5s，时间步长1ms。图6 是图5所示的二维简单介质模型的地震记录剖面；其中，图6(a)第20炮的 观测记录x分量，图6(b)第20炮的观测记录z分量，图6(c)利用传统方 法所得的第20炮的预测记录x分量，图6(d)利用传统方法所得的第20炮 的预测记录z分量，图6(e)利用本发明所得的第20炮的预测记录x分量， 图6(f)利用本发明所得的第20炮的预测记录z分量。由图6所示地震记录 对比可知，利用传统多波角度域共成像点道集所得的预测地震记录与观测地震 记录在同相轴的振幅和相位上有较大的差异，而利用本发明的多波角度域共成 像点道集所得的预测地震记录与观测地震记录一致性较好，间接表明了本发明 所得的多波角度域共成像点道集精度更好。图7是图5所示的二维简单介质模 型的多波局部偏移距域共成像点道集；其中，图7(a)是利用传统方法所得 的纵波局部偏移距域共成像点道
集，图7(b)是利用传统方法所得的横波局 部偏移距域共成像点道集，图7(c)是利用本发明所得的纵波局部偏移距域 共成像点道集，图7(d)是利用本发明所得的横波局部偏移距域共成像点道 集。图8是图5所示的二维简单介质模型的多波角度域共成像点道集；其中， 图8(a)是利用传统方法所得的纵波角度域共成像点道集，图8(b)是利用 传统方法所得的横波角度域共成像点道集，图8(c)是利用本发明所得的纵 波角度域共成像点道集，图8(d)是利用本发明所得的横波角度域共成像点 道集。由图7和图8可知，基于本发明所得的多波角度域共成像点道集具有更 高的分辨率、信噪比及振幅均衡性。上述结果表明了本发明的正确性。
[0162]
实例3：
[0163]
图9是本发明提供的marmousi
‑
2偏移模型；其中，(a)纵波速度模型， (b)横波速度模型，(c)密度模型。该模型是验证各种偏移方法成像效果的 国际标准模型之一。模型深度为5.4km，横向宽度为27.2km。偏移所用空间 网格大小为10m，共计109炮，起始炮点位于模型2.55km处，炮点置于150 m深度处，炮间隔为200m，中间放炮，两边接收，每炮均为501道接收，最 小偏移距0m，最大偏移距2500m，道间距10m，记录时间长度8.0s，时间 步长1ms，采用主频为15hz的雷克子波作为震源时间函数。图10是图9所 示marmousi
‑
2模型的多炮叠加偏移剖面：其中，图10(a)第55炮的观测记 录x分量，图10(b)第55炮的观测记录z分量，图10(c)利用传统方法所 得的第55炮的预测记录x分量，图10(d)利用传统方法所得的第55炮的预 测记录z分量，图10(e)利用本发明所得的第55炮的预测记录x分量，图10(f)利用本发明所得的第55炮的预测记录z分量。由图10所示地震记录 对比可知，利用传统多波角度域共成像点道集所得的预测地震记录与观测地震 记录在同相轴的振幅和相位上有较大的差异，而利用本发明的多波角度域共成 像点道集所得的预测地震记录与观测地震记录一致性较好，间接表明了本发明 所得的多波角度域共成像点道集精度更好。图11是图9所示marmousi
‑
2模型 的纵波局部偏移距域共成像点道集；其中，图11(a)是利用传统方法所得的 纵波局部偏移距域共成像点道集，图11(b)是利用本发明所得的纵波局部偏 移距域共成像点道集。图12是图9所示marmousi
‑
2模型的横波局部偏移距域 共成像点道集；其中，图12(a)是利用传统方法所得的横波局部偏移距域共 成像点道集，图12(b)是利用本发明所得的横波局部偏移距域共成像点道集。 图13是图9所示marmousi
‑
2模型的纵波角度域共成像点道集；其中，图13 (a)是利用传统方法所得的纵波角度域共成像点道集，图13(b)是利用本 发明所得的纵波角度域共成像点道集。图14是图9所示marmousi
‑
2模型的横 波角度域共成像点道集；其中，图14(a)是利用传统方法所得的横波角度域 共成像点道集，图14(b)是利用本发明所得的横波角度域共成像点道集。由 图13和图14可以看出，本发明有提高了角度域共成像点道集的分辨率、信噪 比及振幅随着反射角变化时的能量均衡性。图15是图9所示marmousi
‑
2模型 的角度域共成像点道集叠加偏移剖面：其中，图15(a)是利用传统方法所得 的纵波角度域共成像点道集叠加偏移剖面，图15(b)利用传统方法所得的横 波角度域共成像点道集叠加偏移剖面，图15(c)是利用本发明所得的纵波角 度域共成像点道集叠加偏移剖面，图15(d)是利用本发明所得的横波角度域 共成像点道集叠加偏移剖面。由图15对比可知，本发明所得的多波叠加剖面 质量最高，其反射界面的分辨率高，振幅均衡，断层的断面清晰，断点位置准 确，绕射能量收敛较好，假象较少。上述结果验证了本发明在复杂模型中的有 效性。综上所述，本发明在复杂的地质地球物理模型中具有很好的
可行性及实 用性。
[0164]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是 与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于 实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较 简单，相关之处参见方法部分说明即可。
[0165]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施 例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的 一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变 之处。综上所述，本说明书内容不应理解为对本发明的限制。