一种基于叠前地震射线阻抗反演的不整合圈闭储层岩性预测方法与流程

本发明属于地球物理勘探技术，涉及对不整合圈闭储层描述更可靠的射线阻抗反演和岩性预测方法。

背景技术：

常规的弹性阻抗反演技术适用于与围岩阻抗差异不大的含气砂岩储层的岩性分辨，并在小角度范围内具有较高的精度。虽然弹性阻抗的反演流程较为容易实现，但是精确的叠前地震角道集的构建、混合相位子波估算、高精度反射系数序列反演等关键技术往往对反演结果起到决定性作用。另外，实际数据应用显示，弹性阻抗反演在对非均质性强、储层物性复杂的油气藏进行岩性预测时，常面临诸如数值随角度变化剧烈、大入射角范围内精度不足、预测结果对固有参数(横波纵波速度比)敏感等问题。

与常规的叠前反演技术不同的是，射线波阻抗反演方法及技术能够得到反映岩性参数随频率变化特征的射线波阻抗，可作为岩性和流体识别的预测工具。不整合圈闭储层具有的非均质性强、储层物性复杂等特点，给常规弹性阻抗反演技术带来很多新的困难，叠前射线阻抗反演成为解决该问题的关键技术。

技术实现要素：

本发明的目的是针对现有技术存在的问题,提供一种射线阻抗受入射角范围等影响更小，更适用于描述非均质性强、岩性变化大、储层物性复杂的不整合圈闭储层的基于叠前地震射线阻抗反演的不整合圈闭储层岩性预测方法。

本发明通过如下技术方案实现：

一种基于叠前地震射线阻抗反演的不整合圈闭储层岩性预测方法包括：采用弯曲射线追踪法将经过叠前偏移的共成像点道集从偏移距域转换到射线参数域；对所有的共成像点道集转换后，通过抽取常数射线参数的数据和部分数据叠加构建共射线参数剖面；对共射线参数剖面迭代反演得到高精度地震反射系数序列；在射线参数域中采用广义线性反演算法对反射系数道进行反演求取射线波阻抗；采用能够反映岩性参数随频率变化特征的射线阻抗，结合对其它多种弹性阻抗参数的综合分析，联合探区储层物性特征和地质特征，开展基于统计岩石物理的岩性定量解释方法，实现不整合圈闭储层的岩性预测。

上述方案进一步包括：

(1)将地震采集得到并经过偏移成像的CIP道集通过弯曲射线追踪法，从原本的偏移距域转换到射线参数域。保证了沿着某一条射线参数为常数的射线路径地震反射记录来自于地下同一空间点的反射。根据下面公式通过迭代计算得到共成像点道集中各个反射点的水平方向坐标x_i(i＝0,1,…,2k)

其中u_i为水平慢度，遵从u_i＝u_2k-i对称，其中i＝0,1,…,2k；s_i是点x_i和x_i-1 之间的射线线段长度；

利用sinc方程对射线参数域的波形进行空间插值重构以确保空间连续性和信噪比，得到对应射线参数域的共成像点道集。

(2)对所有的共成像点道集转换后，通过抽取常数射线参数的数据和部分数据叠加构建共射线参数剖面。

(3)通过优化Lp范数和柯西约束两种稀疏反演方法将稀疏约束与地震反射系数分布相结合，通过迭代反演得到高精度地震反射系数序列。

(4)在射线参数域中，将反射系数剖面上的相对振幅值转换成实际的反射系数值后采用广义线性反演算法对反射系数道进行反演求取射线波阻抗。

将从CRP剖面反演得到的反射系数剖面用于反演，求取射线波阻抗。沿着一条具有常数射线参数值的射线路径，其反射系数可以表示为

r(t_i；p)＝f(RI_i(p),RI_i+1(p))， (2)

将该表达式进行Taylor级数展开，并表示为矩阵―向量的形式：

r(RI)≈r(RI₀)+G·(RI-RI₀), (3)

式中RI是待求解的波阻抗剖面，RI₀是波阻抗剖面的初始模型，G是反射系数对波阻抗的偏微分矩阵。求解此一阶Taylor展开式可得到模型修正量δRI＝RI-RI₀。

通过线性迭代方式反演最终得到对应不同射线参数的射线波阻抗。

(5)利用反演得到的对应不同射线参数的对岩性敏感的射线波阻抗，结合其它类型阻抗，基于统计岩石物理技术，我们可以定量地震属性与储层描述之间的关系，合理预测岩性。

本发明所用射线阻抗受入射角范围等影响更小，更适用于描述非均质性强、岩性变化大、储层物性复杂的不整合圈闭储层，且对岩性的预测更加符合地质规律。

附图说明

图1本发明的一种具体实施例的流程图。

图2(a)弯曲射线追踪迭代方法；点线是弯曲射线追踪中的初始路径；实线是该指定同相轴的实际射线路径；虚线是迭代五次后得到射线路径的中间结果。(b)射线参数域共成像点道集。

图3(a)射线参数剖面及对应频谱p＝50ms/km。(b)射线参数剖面及对应频谱p＝100ms/km。(c)射线参数剖面及对应频谱p＝150ms/km。(d)射线参数剖面及对应频谱p＝200ms/km。

图4(a)-(d)分别对应不同射线参数剖面的反射系数反演结果p＝50，100,150和200ms/km。

图5(a)-(d)分别对应射线波阻抗反演结果p＝50，100,150和200ms/km。

图6(a)Monte Carlo模拟估算对应岩性的弹性参数。(b)多变量概率密度函数的建立。(c)某井位声波阻抗与射线阻抗p＝200ms/km交会图及对应不同岩相的二元概率密度函数。(d)不整合圈闭储层岩性预测结果。(e)井位附近沿层二维线岩性预测结果

具体实施方式

本发明通过如下技术方案实现：

(1)将地震采集得到并经过偏移成像的CIP道集，采用在地震层析中普遍使用的方法—弯曲射线追踪法，从原本的偏移距域转换到射 线参数域。保证了沿着某一条射线参数为常数的射线路径地震反射记录来自于地下同一空间点的反射。

跟据费马原理，射线路径S是旅行时T最小的路径。

式中v是速度，s是射线弧长。在每个共成像点道集中，速度场仅沿深度方向变化，待计算的路径可以离散为具有固定深度间隔的折线路径。

S＝{x₀,x₁,…x_k,…x_2k}， (5)

其中每个x代表各个共成像点道集中反射点的水平方向坐标。则旅行时可以表示为

式中u_i＝1/v_i，ds_i是点x_i和x_i-1之间的射线线段长度。在射线追踪中，我们将x^k作为反射点，而射线路径的结束点x₀和x_2k固定。由此费马原理可表示成

它指出了在所有可能的地震射线路径中，地震能量沿其一阶微分为零的路径传播。如果假设偏移后的共成像点道集局部水平，(6)式可微并表示为

其中i＝1,…,2k-1。使用微分项的二阶精度表达式，(9)式给出了一个三对角的线性方程组：

其中共成像点道集中各个反射点的水平方向坐标x_i(i＝0,1,…,2k)，u_i为水平慢度，遵从u_i＝u_2k-i对称，其中i＝0,1,…,2k；s_i是点x_i和x_i-1之间的射线线段长度；

其中每个x代表各个共成像点道集中反射点的水平方向坐标；u_i＝1/v_i；慢度变量遵从u_i＝u_2k-i对称，其中i＝0,…,2k。

跟据以上推导公式通过迭代计算得到未知量x。迭代效果如图2(a)所示。

利用sinc方程对射线参数域的波形进行空间插值重构以确保空间连续性和信噪比

其中A(p)是某一时刻射线参数为p的振幅；x是原始共成像点道集对应的偏移距，x是共成像点道集的空间采样间隔；k表示插值窗的长度，按照重采样定理，k的长度至少为8才能插值后振幅能够保真。得到对应射线参数域的共成像点道集,如图2(b)所示。

(2)对所有的共成像点道集转换后，通过抽取常数射线参数的数据和部分数据叠加构建共射线参数剖面，不同p值的叠前数据具有不同的频带特征，如图3(a)-3(d)所示。具有较小p值的数据具有较高的主频，纵向分辨率更好有利于薄层的分辨。但频带宽度较窄，信噪比 也较低。随着射线参数的增加，剖面的主频逐渐降低，但是带宽逐渐增加。

(3)通过优化柯西约束稀疏反演方法,在稀疏性和反演残差之间采取折衷，修正柯西参数,将稀疏约束与地震反射系数分布相结合，通过迭代反演得到的高精度地震反射系数序列如图4(a)-(d)所示。

(4)在射线参数域中，对每一个CRP剖面推导一个与其p值相对应的标定因子,将反射系数剖面上的相对振幅值转换成实际的反射系数值后采用广义线性反演算法对反射系数道进行反演求取射线波阻抗。

将从CRP剖面反演得到的反射系数剖面用于反演，求取射线波阻抗。沿着一条具有常数射线参数值的射线路径，其反射系数可以表示为

r(t_i；p)＝f(RI_i(p),RI_i+1(p))， (11)

将该表达式进行Taylor级数展开，并表示为矩阵-向量的形式：

r(RI)≈r(RI₀)+G_·(RI-RI₀), (12)

式中RI是待求解的波阻抗剖面，RI₀是波阻抗剖面的初始模型，G是反射系数对波阻抗的偏微分矩阵。求解此一阶Taylor展开式可得到模型修正量δRI＝RI-RI₀。

通过线性迭代方式反演最终得到图5(a)-(d)所示对应不同射线参数的射线波阻抗。

(5)利用反演得到的对应不同射线参数的对岩性敏感的射线波阻抗，结合其它类型阻抗，首先通过蒙特卡洛模拟(图6a)探究不同岩相固有的纵横波速度和密度的变化，钻井未遇地质情况可由岩石物理建模得到，依据二元概率密度函数(图6b)的建立方法计算不同岩相对应的关于地震属性(射线波阻抗与声波阻抗)交会的概率密度函数并经过校正作为分类地震数据的似然函数(图6d)，结合实际地质情况作为约束，根据贝叶斯准则我们可以定量地震属性与储层描述之间 的关系，合理预测岩性(图6d、6e)。