技术领域
[0001] 本
发明涉及地球物理勘探技术，更具体地讲，涉及一种基于P波和Q波属性获取体曲率，最终实现预测裂缝发育的方法和设备。
背景技术
[0002] 据统计，全球石油
天然气产量有一半以上分布于裂缝性储层中。在我国，裂缝性储层也广泛发育，如四川、华北、长庆、塔里木、克拉玛依、胜利、吉林、辽河、青海、玉
门等许多油田都发现了裂缝性油气田。裂缝性储层又以
碳酸盐类为主，然而由于碳酸盐岩储层特有的低孔、非均质性和
各向异性，大大增加了油气勘探的难度。裂缝是此类储层主要油气储集和运移的通道。因此，准确地预测此类储层中的裂缝，就等于直接在此类储层中找到了油气。
[0003] 在常用的裂缝预测方法中，利用
地震信息识别和预测裂缝的方法由于在横向连续性上的优势，因此在裂缝描述和预测中具有不可替代的作用。其
基础是各向异性理论，具体而言，又可分为三大类：
[0004] 一是基于成因分析的预测方法。它主要在从裂缝形成条件出发，如裂缝发育的有利构造部位、有利沉积相带、有利岩性、有利厚度等，间接预测裂缝发育带。
[0005] 二是地震横波预测法。因为横波在裂隙介质中分裂为快、慢横波(横波双折射)，通过对地震资料的坐标变换得到横波的多个分量，这样就可以得到我们所需的信息，进而实现裂缝检测。
[0006] 三是基于裂缝纵波地震响应的预测方法。通过分析裂缝纵波的地震响应特征，如相干体、振幅、波阻抗、吸收系数、
地层曲率、分频属性、属性体和属性差异体，直接预测裂缝的展布特征。
[0007] 基于成因分析的预测方法只能对裂缝进行间接预测；地震横波在裂缝检测上具有相对的优势，但由于成本和处理技术等问题，近期难以推广应用；而纵波地震预测方法则是一种既直接又现实的裂缝预测方法，并且在国内外部分油气田也有成功应用的先例。综合比较而言，基于裂缝纵波地震响应的预测方法更为实际可行。
[0008] 在基于裂缝纵波地震响应的裂缝预测中，包括基于叠前地震资料的裂缝预测和基于叠后地震资料的裂缝预测方法两大类，其中后者又主要包括相干、曲率等方法。
[0009] 目前，曲率计算都是以倾
角和方位角作为输入，且主要采用复数道分析
算法或离散倾角扫描算法。复数道分析算法由于只是用单点的倾角信息，没有平滑掉背景噪声，因此效果较差；而离散倾角扫描算法的效果虽好，但是计算速度太慢。
发明内容
[0010] 本发明目的在于提供体曲率获取方法和设备。
[0011] 为实现上述目的，根据本发明的一方面，提供了一种获取体曲率的方法，包括：A)获取三维地震数据的每个地震数据
采样点的线方向、道方向和时间/深度方向的梯度矢量；B)基于获取的梯度矢量建立地震数据采样点的梯度结构张量；C)对梯度结构张量的元素进行平滑；D)基于平滑后的梯度结构张量获取P波属性和Q波属性，其中，P波属性是一种反应地下岩层沿线方向的变化率的属性，Q波属性是一种反应地下岩层沿道方向的变化率的属性；E)基于P波属性和Q波属性获取体曲率。
[0012] 可选地，步骤C)中可通过下面的高斯
滤波器对梯度结构张量的元素进行平滑：
[0013]
[0014] 其中，xj，yj和zj表示沿着x，y和z轴的待平滑的梯度结构张量的元素为中心的第j道，σ是平滑时窗。
[0015] 可选地，基于平滑后的梯度结构张量获取P波属性和Q波属性可以是对平滑后的梯度结构张量进行特征值分解得到P波属性和Q波属性，其中，P波属性是进行特征值分解后的最大特征值所对应的
特征向量的第一个元素和第二个元素、第三个元素的差的比值；Q波属性是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第二个元素、第三个元素的差的比值。
[0016] 可选地，步骤E)可包括：基于P波属性和Q波属性获取下面的反射面的二次曲面表达式的系数：
[0017] z(x,y)＝ax2+by2+cxy+dx+ey+f
[0018] 其中，a的值是P波属性沿x方向的偏导数值的一半，b的值是Q波属性沿y方向的偏导数值的一半，c是P波属性沿x方向的偏导数值和Q波属性沿y方向的偏导数值的和的一半，d的值等于P波属性的值，e的值等于Q波属性的值，f是常数；基于获取的二次曲面表达式的系数获取体曲率。
[0019] 可选地，所述方法可还包括：读取三维地震数据，其中，三维地震数据是叠后地震数据。
[0020] 根据本发明的另一方面，提供了一种获取体曲率的设备，包括：梯度矢量获取模
块，用于获取三维地震数据的每个地震数据采样点的线方向、道方向和时间/深度方向的梯度矢量；梯度结构张量获取模块，用于基于获取的梯度矢量建立地震数据采样点的梯度结构张量；平滑模块，用于对梯度结构张量的元素进行平滑；P波属性和Q波属性获取模块，用于基于平滑后的梯度结构张量获取P波属性和Q波属性，其中，P波属性是一种反应地下岩层沿线方向的变化率的属性，Q波属性是一种反应地下岩层沿道方向的变化率的属性；体曲率获取模块，用于基于P波属性和Q波属性获取体曲率。
[0021] 可选地，平滑模块可通过下面的高斯滤波器对梯度结构张量的元素进行平滑：
[0022]
[0023] 其中，xj，yj和zj表示沿着x，y和z轴的待平滑的梯度结构张量的元素为中心的第j道，σ是平滑时窗。
[0024] 可选地，基于平滑后的梯度结构张量获取P波属性和Q波属性可以是对平滑后的梯度结构张量进行特征值分解得到P波属性和Q波属性，其中，P波属性是进行特征值分解后的最大特征值所对应的特征向量的第一个元素和第二个元素、第三个元素的差的比值；Q波属性是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第二个元素、第三个元素的差的比值。
[0025] 可选地，体曲率获取模块可包括：系数获取单元，用于基于P波属性和Q波属性获取下面的反射面的二次曲面表达式的系数：
[0026] z(x,y)＝ax2+by2+cxy+dx+ey+f
[0027] 其中，a的值是P波属性沿x方向的偏导数值的一半，b的值是Q波属性沿y方向的偏导数值的一半，c是P波属性沿x方向的偏导数值和Q波属性沿y方向的偏导数值的和的一半，d的值等于P波属性的值，e的值等于Q波属性的值，f是常数；体曲率计算单元，用于基于获取的二次曲面表达式的系数获取体曲率。
[0028] 可选地，所述设备可还包括：地震数据读取模块，用于读取三维地震数据，其中，三维地震数据是叠后地震数据。
[0029] 通过本发明，将实现一种基于P波和Q波属性的多尺度体曲率裂缝预测方法，在具有较高的计算性能的同时，能清晰地刻画裂缝发育细节，抗噪能
力更强，帮助识别油气储集和运移的通道，从而精确地确定勘探目标和降低勘探
风险，对裂缝性储藏的勘探起到重要的作用。
[0030] 将在接下来的描述中部分阐述本发明另外的方面和/或优点，还有一部分通过描述将是清楚的，或者可以经过本发明的实施而得知。
附图说明
[0031] 通过下面结合附图进行的详细描述，本发明的上述和其它目的、特点和优点将会变得更加清楚，其中：
[0032] 图1示出了根据本发明示例性
实施例的获得体曲率的方法的
流程图；
[0033] 图2示出了根据本发明示例性实施例的获得体曲率的设备的
框图。
具体实施方式
[0034] 这里公开详细示出的示例实施例。然而，这里公开的特定结构和功能细节仅是有代表性的，目的是描述示例实施例。然而，可以以多种
变形实施本发明，并且本发明不能被理解为仅限于在此阐述的示例实施例。
[0035] 图1示出了获得体曲率的方法的流程图。
[0036] 参照图1，在步骤S101，获取三维地震数据的每个地震数据采样点的线方向、道方向和时间/深度方向的梯度矢量。
[0037] 三维地震数据是叠后地震数据，例如，以本领域公知的方法进行野外勘探采集，可得到原始采集的数据，然后，经过地震资料处理(包括本领域公知的球面发散校正、对地层吸收的Q补偿、振幅处理、子波反褶积、地表一致性静校正、速度分析、动校正、剩余静校正、
水平
叠加、偏移等步骤)可得到叠后地震数据。
[0038] 例如，梯度矢量可包括线方向、道方向和时间/深度方向的三个梯度分量。可通过求取沿着线方向、道方向和时间/深度方向的偏导数的方式获得梯度矢量。可分别通过下面的式(1)计算地震数据三个方向 和 的梯度，从而得到梯度矢量g，
[0039]
[0040] 其中，x、y和z分别是对应于线方向、道方向和时间/深度方向的坐标，u是沿着对应方向的地震采样数据离散序列。
[0041] 在步骤S102，基于获取的梯度矢量建立地震数据采样点的梯度结构张量。
[0042] 作为一个示例，可对式(1)中的梯度向量g做外积的方式(即梯度向量和其转置相乘)来建立式(2)所示的地震数据采样点的梯度结构张量T，梯度结构张量的表达式如下：
[0043]
[0044] 其中，g1、g2和g3是梯度矢量的3个元素。
[0045] 在步骤S103，对梯度结构张量的元素进行平滑。
[0046] 作为一个示例，可采用高斯核平滑函数分别对张量的每个元素进行平滑。
[0047] 高斯滤波器的表达式为：
[0048]
[0049] 其中，xj，yj和zj分别表示沿着x，y和z轴的以待求导的分析点为中心的第j道，σ是平滑时窗。当地震数据经过最佳采样时，那么三个方向的理想尺度因子σ分别可设为时间轴采样(sample)间隔、
共中心点集(CMP)间隔和线间距(line spacing)大小的三倍。
[0050] 平滑后的梯度结构张量TS如式(3)所示：
[0051]
[0052] 其中，Ts是平滑后的梯度结构张量，G是高斯滤波器，g1、g2和g3是梯度矢量的3个元素。
[0053] 在步骤S104，基于平滑后的梯度结构张量获取P波属性和Q波属性。其中，P波属性是一种反应地下岩层沿线方向的变化率的属性，Q波属性是一种反应地下岩层沿道方向的变化率的属性。
[0054] 作为一个示例，可对平滑后的梯度结构张量TS做如下的特征值分解：
[0055] TSv＝Λv
[0056] 其中，Λ表示对角矩阵，对角线上的值为特征值，分别为λ1，λ2，λ3，它们是由大到小排列的。v表示特征向量矩阵，列向量表示以上特征值对应的特征向量，特征值按照从大到小排列分别为v1，v2，v3。
[0057] 作为一个示例，P波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第一个元素和第三个元素的比值；Q波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第三个元素的比值。例如，可定义P波属性为：
[0058] p＝v1(1)/v1(3)，
[0059] 可定义Q波属性为：
[0060] q＝v1(2)/v1(3)。
[0061] 作为另一个示例，P波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第一个元素和第二个元素、第三个元素的差的比值；Q波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第二个元素、第三个元素的差的比值。可定义P波属性为：
[0062] V1(1)/(V1(2)-V1(3))
[0063] 可定义Q波属性为：
[0064] V1(2)/(V1(2)-V1(3))
[0065] 其中，关于参数的含义，以下面示例进行说明，V1表示进行特征值分解后的最大的特征值所对应的特征向量，V3表示进行特征值分解后的最小的特征值所对应的特征向量，V2表示进行特征值分解后的介于最大特征值和最小特征值之间的特征值所对应的特征向量，所述的每个特征向量包括3个元素，V1(2)表示最大的特征值所对应的特征向量的第2个元素。
[0066] 在步骤S105，可基于P波属性和Q波属性获取体曲率。
[0067] 作为一个示例，可采用
现有技术中常用的二次曲面，例如，可将反射面的二次曲面表述为：
[0068] z(x,y)＝ax2+by2+cxy+dx+ey+f
[0069] 那么，基于P波属性和Q波属性和如下表达式获取二次曲面的系数，Dx(Dxz)＝2a＝Dxp；Dy(Dyz)＝2b＝Dyq；Dxp+Dyq+＝2c；d＝p；e＝q。a的值是P波属性沿x方向的偏导数值的一半，b的值是Q波属性沿y方向的偏导数值的一半，c是P波属性沿x方向的偏导数值和Q波属性沿y方向的偏导数值的和的一半，d的值等于P波属性的值，e的值等于Q波属性的值，f是常数。
[0070] 基于现有技术中的体曲率计算公式来获得体曲率，体曲率计算公式表征如下：
[0071] 1)平均曲率
[0072]
[0073] 2)高斯曲率
[0074]
[0075] 3)极大曲率
[0076]
[0077] 4)极小曲率
[0078]
[0079] 5)最大正曲率
[0080]
[0081] 6)最小负曲率
[0082]
[0083] 7)形态指数
[0084]
[0085] 8)倾向曲率
[0086]
[0087] 9)走向曲率
[0088]
[0089] 10)等值线曲率
[0090]
[0091] 11)弯曲度
[0092]
[0093] 其中，a、b、c、d、e分别是二次曲面表达式的系数，根据以上曲率计算公式，分别基于a、b、c、d、e计算即可得各体曲率。
[0094] 图2示出了获得体曲率的设备的框图。
[0095] 参照图2，获得体曲率的设备包括：梯度矢量获取模块201、梯度结构张量获取模块202、平滑模块203、P波属性和Q波属性获取模块204和体曲率获取模块205。
[0096] 梯度矢量获取模块201用于获取三维地震数据的每个地震数据采样点的线方向、道方向和时间/深度方向的梯度矢量。具体实施方式和步骤101所述内容一致，在此不再赘述。
[0097] 梯度结构张量获取模块202可用于基于获取的梯度矢量建立地震数据采样点的梯度结构张量。具体实施方式和步骤102所述内容一致，在此不再赘述。
[0098] 平滑模块203可用于对梯度结构张量的元素进行平滑。平滑模块可通过下面的高斯滤波器对梯度结构张量的元素进行平滑：
[0099]
[0100] 其中，xj，yj和zj表示沿着x，y和z轴的待平滑的梯度结构张量的元素为中心的第j道，σ是平滑时窗。
[0101] P波属性和Q波属性获取模块204可用于基于平滑后的梯度结构张量获取P波属性和Q波属性，其中，P波属性可以是一种反应地下岩层沿线方向的变化率的属性，Q波属性可以是一种反应地下岩层沿道方向的变化率的属性。
[0102] 作为一个示例，P波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第一个元素和第三个元素的比值；Q波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第三个元素的比值。
[0103] 作为另一个示例，P波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第一个元素和第二个元素、第三个元素的差的比值；Q波属性可以是进行特征值分解后的最大特征值所对应的特征向量的第二个元素和第二个元素、第三个元素的差的比值。
[0104] 体曲率获取模块205可用于基于P波属性和Q波属性获取体曲率。可基于P波属性和Q波属性获取下面的反射面的二次曲面表达式的系数：
[0105] z(x,y)＝ax2+by2+cxy+dx+ey+f
[0106] 其中，a的值是P波属性沿x方向的偏导数值的一半，b的值是Q波属性沿y方向的偏导数值的一半，c是P波属性沿x方向的偏导数值和Q波属性沿y方向的偏导数值的和的一半，d的值等于P波属性的值，e的值等于Q波属性的值，f是常数；可基于获取的二次曲面表达式的系数获取体曲率。
[0107] 此外，所述设备还可包括地震数据读取模块(未示出)，地震数据读取模块可用于读取三维地震数据，其中，三维地震数据是叠后地震数据。例如，以本领域公知的方法进行野外勘探采集，可得到原始采集的数据，然后，可经过地震资料处理(包括本领域公知的球面发散校正、对地层吸收的Q补偿、振幅处理、子波反褶积、地表一致性静校正、速度分析、动校正、剩余静校正、水平叠加、偏移等步骤)得到叠后地震数据。
[0108] 尽管已经参照其示例性实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离
权利要求所限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。
