高分辨率转换波裂缝预测方法与流程

本发明属于转换波地震资料处理技术领域，具体来讲，涉及一种高分辨率转换波裂缝预测方法。

背景技术：

目前，相干体裂缝特征描述是探寻裂缝性油气藏的一项重要手段，而获取裂缝预测所需高分辨率数据体、裂缝参数的求取一直是裂缝油气藏勘探中重点研究的问题，尤其转换波低信噪比和分辨率的特点，需要提高其分辨率，便于后续解释的薄层分辨。但是，最小相位和白噪假设为前提的常规高分辨率处理方法，在主频带范围限制下，提高转换波分辨率的能力有限，无法达到后续裂缝预测分辨多薄层最小厚度的要求，所以需要采取新的高分辨率处理思路，即：转换波高精度反射系数求取与高频子波高分辨率重构的新思路，该思路不受任何假设和频带限制，能完成常规手段无法分辨多薄层最小厚度的要求，获取裂缝预测满意的高品质剖面。

通常，反射系数的求取需要构建多薄层目标函数，然而，由于常规方法构建的目标函数存在诸多不合理之处，故而，在求取过程中的反射系数存在极值和位置、大小上的误差，这种误差导致后续高分辨率剖面的错像。此外，目标函数的求解方法有很多，比如模拟退火法、共轭梯度法、匹配追踪法等这些算法，虽然都能求取目标函数的反射系数，但精度和稳定性较差，且对新构建的目标函数适应性较差。

三维相干体技术是20世纪90年代后期兴起的一项十分有效的地震解释技术，该技术主要从相邻地震道相互之间的相干性出发，给出一种定量描述，它能突出断层处地震波的变化情况，形成相干值异常区域，从而预测断裂展布。M.Bahorich等提出在振幅上进行互相关相干分析算法，但该算法计算出的断点非相干性点模糊，噪声重、能量弱。K.J.Marfurt年提出基于相似系数相干算法具有稳定性好、断层刻画精度高的特点，但横向分辨率低的劣势无法改变。Gersztenkorn等提出基于本征结构分析的算法，利用Hilbert变换计算相干体提高其稳定性和抗噪性，但该算法用于高陡复杂地区时，如不考虑倾角则不如第二代算法效果好，如考虑倾角则计算效率无法应用于生产中。

技术实现要素：

本发明的目的在于解决现有技术存在的上述不足中的至少一项。例如，本发明的目的之一在于提供一种具有更高分辨率的转换波裂缝预测方法，以使得提取的裂缝预测剖面分辨率更高、裂缝走向更加清晰、细节刻画更加明显。

为了实现上述目的，本发明提供一种高分辨率转换波裂缝预测方法。所述方法包括：采集转换波地震数据，进行预处理和叠前偏移后得到转换波偏移剖面，在转换波偏移剖面上取n道地震记录x(x,t)，x＝1,2,…,n；在偏移剖面上建立稀疏约束和波阻抗约束的多薄层目标函数；对多薄层目标函数利用频域映射与L1范数联合优化算法求取高精度反射系数；构建高频子波同时保持该高频子波具有确定的频率宽度，用高频子波与高精度反射系数得到高分辨率转换波偏移剖面；利用改进的第三代特征值相干体计算法，对高分辨率转换波偏移剖面进行相干数据体切片提取，从而实现裂缝预测；其中，

所述目标函数为下式(1)：

其中，r表示反射系数，a_e代表反射系数的偶分量比例，a_o代表反射系数的奇分量比例，r_e为反射系数的偶分量，r_o为反射系数的奇分量，A为变换矩阵，b为地震数据和子波相关矩阵，λ为稀疏因子，为阻抗因子，C为积分算子矩阵，ξ_p为纵波波阻抗矩阵。

与现有技术相比，本发明的有益效果包括：能够基于转换波地震数据构建更加合理的多薄层目标函数，并能有效提高目标函数求解精度和抗噪性，能够获得高精度反射系数，能够使得提取的裂缝预测剖面分辨率更高、裂缝走向更加清晰、细节刻画更加明显。

附图说明

图1示出了通过现有的预处理和叠前偏移后得到的转换波偏移剖面。

图2示出了根据本发明的方法求取高精度反射系数对应的剖面。

图3示出了根据本发明的方法得到的高分辨率转换波偏移剖面。

图4A示出了图1的转换波偏移剖面的振幅谱数据；图4B示出了根据本发明的方法得到的高分辨率转换波偏移剖面的振幅谱数据。

图5A示出了利用改进的第三代特征值相干体计算法对图1的转换波偏移剖面进行相干切片提取的效果图；图5B示出了根据本发明的方法对高分辨率转换波偏移剖面进行相干切片提取的效果图。

具体实施方式

在下文中，将结合附图和示例性实施例来详细说明本发明的高分辨率转换波裂缝预测方法。

在本发明的一个示例性实施例中，高分辨率转换波裂缝预测方法可通过下面的步骤(1)至步骤(5)来实现：

(1)对采集到的转换波地震数据，经过常规预处理和叠前偏移后得到的转换波偏移剖面，取n道地震记录x(x,t)，x＝1,2,…,n；其中，x表示道数，t表示时间。

(2)在偏移剖面上建立由稀疏约束和波阻抗约束双重约束的多薄层目标函数；

步骤(2)可按下列方式完成：

(a)对于不含噪声的褶积模型，地震合成记录在频率域中可表示为：

S(f)＝W(f)·R(f) (1)

其中，S(f)表示频率域下地震记录(即，上面的x(x，t))，W(f)表示频率域下地震子波，R(f)表示频率域下反射系数；f为频率，时间t经过傅里叶变换后变换到频率f中。

利用频率域下地震记录与地震子波的丰富信息，利用合适的常规反演算法(如，共轭梯度法、匹配追踪法、模拟退火等)消除子波(即，地震子波)影响(如，子波主频、旁瓣干扰等会影响反演求取反射系数精度的内容)，从记录中得到所需的反射系数(这里，为了得到反射系数R(f))。利用反射系数奇偶分解原理，可得到频率域下目标函数表达式：

其中，a_e代表反射系数R(f)的偶分量比例，a_o代表反射系数R(f)的奇分量比例；Re表示R(f)的实部，Im表示R(f)的虚部。

(b)为建立更符合实际地层情况的反射系数模型，对步骤(a)中的多薄层反射系数模型下的目标函数进行推导或改进。

由奇偶反射系数分解原理，得多薄层反射系数奇偶分量可表示为：

r_e(t)代表反射系数r(t)的偶分量，r_o(t)代表反射系数r(t)的奇分量，t表示时间，T表示时间厚度；r(t)为时间域反射系数，r(t)经过傅里叶变换后到上面的频率域反射系数R(f)；g(t)当前分析位置反射系数，g_e(t)代表g(t)的偶分量，g_o(t)代表g(t)的奇分量。

对式(3)进行傅里叶变换，得到其频率域实部与虚部表达式：

其中，π为圆周率。

利用式(2)和式(4)得到多薄层反射模型的目标函数表达式为：

其中，r_e代表反射系数r(t)的偶分量，r_o代表反射系数r(t)的奇分量，Δt表示时移量；i无物理意义，是指数与正弦函数变化标识符。

为了便于求解，将式(5)可写成矩阵形式：

其中，

其中，t为时间，t₁、t₂、……、t_N分别表示第1、2、……、N个时间，f为频率，f₁、f₂、……、f_M分别表示第1、2、……、M个分析频率范围内的频率，T₁、T₂……、T_N/2分别表示第1、2、……、N/2个时间厚度，Δt表示时移量

(c)由于该目标函数(6)(即，式(6))有多个全局最优解，为了降低求解问题时的多解性和不确定性，需要加上稀疏化约束，将目标函数改为范数求解形式为：

其中，λ为稀疏因子，λ一般在0-1范围内取值，|| ||₁表示L1范数，表示欧几里得(Euclidean)范数的平方。

(d)由于波阻抗可看作是反射系数对时间的积分，因此发明人把它引入到目标函数中作为反射系数r的先验约束。

由于转换波相对于纵波来说其信噪比和分辨率较低,为了提高反射系数精度,本发明使用分辨率更高的纵波波阻抗作为转换波的约束条件。

若I(t)为纵波波阻抗值，则纵波相对波阻抗ε_t为：

其中，I(t₀)表示纵波初始波阻抗值，r为反射系数，η表示在时间t₀与时间t中取值符合，t₀表示初始时间，t_M-1表示结束时间。

把上式(8)写成矩阵形式：

ξ_p＝Cr (9)

其中C为积分算子矩阵，其离散形式可表示为：

则由最小平方定义波阻抗约束得到新的目标函数为：

其中，r表示反射系数，a_e和a_o分别代表反射系数r(t)的偶奇分量比例，r_e是反射系数的偶分量，r_o是反射系数的奇分量；A为变换矩阵，b是地震数据和子波相关矩阵，λ是稀疏因子，是阻抗因子(其大于零即可)，C为积分算子矩阵，ξ_p是纵波波阻抗矩阵)

(3)对于步骤(2)中的式(11)，利用频域映射与L1范数联合优化算法求取高精度反射系数。

本步骤(3)可按下列方式完成：

(a)对式(11)进行求解中，由于波阻抗约束项是由测井信号得到的已知先验约束条件，所以式(11)的求解类似于求矩阵方程Ax＝b在含有噪声情况下的稀疏解。

可将目标函数(11)的求解问题转化为下式：

其中，υ为极小权重值，G为复数集合；本式(12)中，A表示A_e时b表示b_e，X表示r_e；A表示A_o时b表示b_o，X表示r_o。

如果引入变量q∈G^m，那么式(12)等价于式(13)：

式(13)对应的增广拉格朗日子问题可表示为：

其中，y^T表示拉格朗日乘子的共轭转置运算，β为罚参数。

利用式(14)进行频域映射与L1范数联合优化算法实现步骤如下：

1)令k＝0对q^k，x^k，y^k赋初始值，并给定υ，β常数值(这里，υ，β为人为给定，可为大于零的常值)，然后进行以下步骤2)至5)的算法运算，如果满足终止准则(例如，终止准则可以为：最优解值或者迭代次数)，则算法完成运算；否则，进行步骤2)；

2)令x＝x^k，y＝y^k，求解r的子问题得到q^k+1：

3)令q＝q^k+1，y＝y^k则关于x的极小化问题式(14)等价于式(16)：

那么可以通过式(17)近似求解来完成对式(16)的精确求解：

其中，τ为大于零的邻近参数(这里，τ为人工给定，其大于零即可)，g^k为x＝x^k时二次项的梯度，那么式(17)通过下式(18)求解得到x^k+1：

其中，ο表示逐元素相乘，Shrink(,)表示一维收缩算子，sign表示符号函数。

4)令x＝x^k+1，q＝q^k+1，求解y的子问题得到y^k+1：

y^k+1＝y^k-γβ(Ax^k+1`+q^k+1-b) (19)

其中，γ为大于零的常数(这里，γ人工给定，其大于零即可)。

5)令k＝k+1，重复步骤1)～步骤4)进行迭代运算。

(4)基于步骤(3)反射系数与高频子波重构，得到高分辨率转换波偏移剖面。

构建高频子波的同时需要保持该高频子波具有一定的频带宽度，以使得高分辨率转换波偏移剖面是宽频的。可利用常规的精度较高的俞式子波建立高频子波信号，亦可用诸如雷克子波等零相位子波建立高频子波信号。

时间域中俞氏子波表达式为：

其中q、p分别表示子波的低频值和高频值，t为时间。

用式(20)的子波与步骤(3)求取的高精度反射系数得到高分辨率转换波偏移剖面。

(5)对步骤(4)的高分辨率转换波偏移剖面进行改进的第三代特征值相干算法C3相干体切片裂缝预测。

相干数据体技术主要利用相邻地震道的相似性原理，描述地层及岩性的横向不均匀性。通过相干数据体算法，对地震数据体的不连续性进行系统分析，使解释人员能有效识别构造和断层的分布情况，避免解释的随意性。

(a)为了尽可能压制噪声干扰，提高相干体计算数据的分辨率，Gersztenkorn等提出了第三代特征值相干体计算法。

首先，利用Marfurt等在第二代多道相似性相干体计算中定义的椭圆和矩形分析窗口中，以坐标为(x,y)的分析点为中心，时窗中心t＝nΔt，则时间窗口内含2k+1个元素的分析道与窗口内的J道数据排列成矩阵：

其中，表示第j道在时间t＝mΔt-pΔx_j-qy_j沿着视倾角(p,q)插值后的振幅值。

那么矩阵U的第m列是分析体内每个地震道在时间t＝mΔt-pΔx_j-qy_j(j＝1,2,…,J)的振幅值，则沿着一对视倾角(p,q)，时窗中心在t＝mΔt的分析窗口内的2k+1个采样点的数据协方差矩阵为：

如λ_j为分析窗口内地震数据协方差矩阵式(22)的特征值，且有λ₁≥λ₂≥…≥λ_J。

则在以分析点为中心的分析窗口内，沿着视倾角(p,q)的J道c₃特征值相干计算算法式可表示为：

则在以分析点为中心的分析窗口内，沿视倾角(p,q＝0)的J道c₃特征值相干计算算法式可表示为：

则沿着所有倾角的c₃特征值相干计算算法式可表示为：

(b)因为第三代特征值相干算法C3是在子空间中计算，比前两代相干算法更能很好地压制噪声，但其计算量大、效率慢的缺点却影响了其在实际生产中的经济效益。

对该算法进行常规改进。例如，可利用基于相似性算法的高效率性估算各点的视倾角和方位角，然后在大于10倍的分析窗口内计算视倾角和方位角的平均值再在以分析点为中心的分析窗口内，沿着区域平滑后的视倾角计算地震数据协方差矩阵的特征值及相干值：

利用式(26)对步骤(4)的高分辨率剖面进行相干数据体切片提取，从而进行裂缝预测。

图1示出了通过现有的预处理和叠前偏移后得到的转换波偏移剖面，其横坐标表示CMP道号，图1中也示出了测井系数(即，图1中位于道号791和841之间的纵向曲线)。图2示出了根据本发明的方法求取高精度反射系数对应的剖面。由图2可以看出，本发明的方法求取高精度反射系数与测井系数(纵向曲线)匹配较好，各层位置和极性求取较佳。

图3示出了根据本发明的方法得到的高分辨率转换波偏移剖面。可以看出，图3明显比图1的分辨率高，且浅、中、深层的同相轴(例如，800ms附近的浅层、1000ms附近的中层、1150ms附近的深层的同相轴)连续性更好、断点更清晰、断层更明了，便于后续解释人员的裂缝预测处理。

图4A示出了图1的转换波偏移剖面的振幅谱数据；图4B示出了根据本发明的方法得到的高分辨率转换波偏移剖面的振幅谱数据。由图4A和图4B可知，经本发明的方法处理而得到的高分辨率转换波偏移剖面的振幅谱数据的频带得到明显拓宽，而且其振幅谱低频信息得到完整的保留。

图5A示出了利用改进的第三代特征值相干体计算法对图1的转换波偏移剖面进行相干切片提取的效果图；图5B示出了根据本发明的方法对高分辨率转换波偏移剖面进行相干切片提取的效果图。由图5A和图5B可知，经本发明的方法处理后的切片相干效果更好，裂缝刻画细节更清晰，避免了低分辨率数据为后续解释人员带来的假象。

综上所述，本发明通过利用稀疏约束和波阻抗双重约束构建多薄层目标函数，提高目标函数求解精度和抗噪性；利用频域映射与L1范数联合优化算法求取高精度反射系数，并通过与高频子波重构得到分辨率更高的新地震剖面，改善转换波的分辨率效果，利用改进的C3相干体切片进行裂缝预测，可有效提高地震勘探裂预测和描述地下构造的能力。也就是说，本发明的高分辨率转换波裂缝预测方法能够使得提取的裂缝预测剖面分辨率更高、裂缝走向更加清晰、细节刻画更加明显。

尽管上面已经结合附图和示例性实施例描述了本发明，但是本领域普通技术人员应该清楚，在不脱离权利要求的精神和范围的情况下，可以对上述实施例进行各种修改。