一种基于等级相关分析的多分量转换波裂缝预测方法

摘要

本发明提供了一种基于等级相关分析的多分量转换波裂缝预测方法，所述预测方法包括以下步骤：选取计算时窗，将待处理采样点的径向分量和横向分量按照计算时窗进行截取；对待处理采样点进行等级相关分析，得到最大等级相关系数；判断待处理采样点处是否具有裂缝发育，确定裂缝发育方位；对判定有裂缝发育的待处理采样点径向分量和横向分量进行分离，得到快横波和慢横波；根据快横波和慢横波时差计算各向异性系数；重复所述选取计算时窗步骤至所述计算各向异性系数步骤，循环至下一个待处理采样点直到地震数据计算完毕。根据本发明的方法采用了判别方法对计算结果进行判别，剔除裂缝明显不发育区域，在对裂缝发育进行统计时更为准确。

说明书

技术领域

本发明涉及石油勘探领域，具体来讲，涉及一种基于等级相关分析的多分量转换波裂缝预测方法。

背景技术

裂缝性油气储层是当今石油地球物理勘探研究的重点之一，如何精确地判定裂缝在地下的存在状态又是重点中的难点。从世界上已开发油气田统计数字来看，裂缝型储集层在油气资源和生产能力方面，约占世界总量的一半。在油藏评价、区块的评估和提高油气采收率中，识别和刻画断层与裂缝系统是至关重要的步骤。在当前裂缝预测方法中，常用的多为利用叠后地震数据计算相干或曲率等属性或属性体预测裂缝，它们多适用于预测大尺度断裂带，无法充分利用叠前信息及地层中的各向异性信息，预测结果不够精细。因此，需要探索更为精细的裂缝预测技术。

随着多波研究的不断深入，多波技术越来越多地应用到油气勘探开发中。横波分裂作为天然地震、石油勘探领域的一项具有强大潜力的技术被地质学家和地球物理学家广泛关注。横波分裂理论即当横波穿过方位各向异性介质时会形成双折射，分裂成沿平行于裂隙方向传播的S1波和沿垂直裂隙方向传播的S2波，它们各自有对应的速度，前者传播的速度较快，叫快横波，后者传播的速度较慢，叫慢横波。由于传播的速度不同，在各向异性介质中形成快慢横波时差、振幅衰减等，这些特性都与各向异性性质有关。因此可根据快、慢波的时差、波形、振幅衰减等反过来研究裂缝系统的方位和密度。这个理论对预测与油气密切相关的裂隙发育带有着十分重要的作用，因此，利用横波分裂预测储层裂缝成为一种直接可靠的方法。

横波分裂的方位和时差是研究方位各向异性和预测裂缝发育带的主要两个参数，因此，准确地计算出裂缝发育方位和时差对于准确刻画裂缝系统在地下的分布形态具有重要意义。目前，常用的方法有互相关法和能量比法。互相关法是Naville(1986)提出的一套根据互相关原理确定横波分裂方位的方法，通过对不同分量数据旋转后求解互相关系数来确定裂缝方位。能量比法是Granger(2002)提出的一种方法，该方法认为当径向分量和横向分量旋转到裂缝发育的方位时，两者的能量比值最大。能量比检测方法当方位角间隔较大时，对裂缝走向的检测精度不足，尤其对于浅层的EDA介质，由于不规则、不均匀的分方位角道集分布，计算误差较大，因此只适合宽方位角三维观测系统。以上两种方法均对计算结果没有一个过滤机制，即一些裂缝不发育的地区也计算出了方位和时差，并绘制在裂缝发育图中。并且，由于横波震源比较昂贵加之横波震源激发的横波信噪比和分辨率都较低，利用纯横波震源来进行横波分裂勘探还难以实现大规模应用。然而，多分量转换波勘探克服了纯横波勘探激发难、成本高、静校正量大等缺陷，同时具有信息丰富、兼有纵波和横波的优势、资料信噪比较高等优点，使得多分量转换波勘探成为油气储层探测的有力工具。

发明内容

针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的一个或多个问题。

为了实现上述目的，本发明提供了一种基于等级相关分析的多分量转换波裂缝预测方法，所述预测方法可以包括以下步骤：选取一个包含目标深度范围的计算时窗，将待处理采样点的径向分量地震数据和横向分量地震数据按照计算时窗进行截取；对待处理采样点进行等级相关分析，得到等级相关系数最大值；判断待处理采样点处是否具有裂缝发育，若有裂缝发育，确定裂缝发育方位；对判定有裂缝发育的待处理采样点进行快横波和慢横波分离；根据快横波和慢横波时差计算各向异性系数；重复所述选取计算时窗的步骤至所述计算各向异性系数的步骤，循环至下一个待处理采样点直到地震数据计算完毕，其中，所述对待处理采样点进行等级相关分析的步骤可包括：将截取的待处理采样点的径向分量地震数据和横向分量地震数据分别记录到数组R和数组T中；对数组R和数组T进行旋转角度θ_j和时移τ_j的变换，得到变换后的分量LR(θ_j,t)和LT(θ_j,t+τ_j)，其中，旋转角度变换式为：

LT(θ_j,t+τ_j)＝LT(θ_j,t)*τ_j，其中，θ_j＝θ_k+Δθ，τ_j＝τ_k+Δτ，θ_k为起始旋转角度，τ_k为起始时移，Δθ为旋转角度循环步长，Δτ为时移循环步长，t为时间；

分别求取LR(θ_j,t)分量和LT(θ_j,t+τ_j)分量的秩，记为ZR_i和ZT_i；利用式1计算等级相关系数，得到等级相关系数最大值r_max＝max{r(θ_j,τ_j)}，其中，式1为其中，D_i＝ZR_i-ZT_i，n为分量中数据个数。

与现有技术相比，本发明的方法采用了判别方法对计算结果进行判别，剔除裂缝明显不发育区域，剔除冗余的背景信息，在对裂缝发育进行统计时更为准确，进行裂缝发育图绘制时使裂缝带的分布和发育情况更加突出。同时，本发明方法提出了一套优化的多分量转换波裂缝预测流程，对裂缝预测过程中的多个步骤进行了方法的优选，包括将等级相关理论引入到多分量转换波裂缝预测中，替代传统的互相关计算，规避了互相关对数据分布的要求，提高了计算的准确度；采用分析-判别模式对裂缝预测结果做进一步精细判断；采用多种模式进行裂缝发育的综合绘制以及井周裂缝发育的统计功能。

附图说明

通过下面结合附图进行的描述，本发明的上述和其他目的和特点将会变得更加清楚，其中：

图1示出了根据本发明示例性实施例的基于等级相关分析的多分量转换波裂缝预测方法流程图。

图2示出了根据本发明示例性实施例的坐标旋转示意图。

图3示出了根据本发明示例性实施例的横波分裂与裂缝发育方位关系示意图。

具体实施方式

在下文中，将结合附图和示例性实施例详细地描述根据本发明的基于等级相关分析的多分量转换波裂缝预测方法。

具体来讲，转换横波是指野外勘探时采用纵波激发，地震波倾斜入射到弹性分界面时会产生反射横波，其中的上行横波为转换横波。转换横波拥有横波的全部性质，当上行转换横波与裂缝走向斜交时，同样会发生分裂现象，分裂为快横波和慢横波。分裂后的快、慢波传播到地表，按照地表观测系统坐标系分解并重新合成，被多分量检波器接收。因此，当测线与裂缝方向不平行时，地面X、Y分量记录中都含有快、慢波信息，然后利用多分量记录进行快、慢波分离，根据快波与慢波的时差、波形、振幅衰减等反过来研究裂缝系统的方位、密度和时差等信息，并根据以上信息预测裂缝发育。

根据横波分裂理论，当横波穿过方位各向异性介质时会形成双折射，分裂成快横波和慢横波。在本发明的裂缝预测的方法中，方位各向异性介质可以是扩容各向异性介质(EDA)，扩容各向异性介质可近似表示由于构造应力产生空间排列垂直裂缝群体引起的各向异性。EDA介质是一种典型的方位各向异性介质，由平行的垂直裂隙、微裂隙或定向的孔隙所引起。EDA介质具有以下特点：(1)裂缝走向与地层最小应力方向垂直，与最大应力方向平行；(2)在地下足够深处，由于地层压力，最小应力通常是水平的，所以裂缝往往是垂直的；(3)偏振波沿应力大的方向传播的速度比沿应力小的方向传播速度大。

图1示出了根据本发明示例性实施例的基于等级相关分析的多分量转换波裂缝预测方法流程图。图2示出了根据本发明示例性实施例的坐标旋转示意图。图3示出了根据本发明示例性实施例的横波分裂与裂缝发育方位关系示意图。

利用本发明的方法进行裂缝预测时，需要用到三种坐标系，分别是野外采集坐标、径向-横向坐标和自然坐标。

(1)野外采集坐标/X-Y坐标系

X轴沿Inline(主测线)方向，以地震排列站号增加方向为正；Y轴沿Crossline方向；Z方向垂直向下。

(2)径向-横向坐标系/R-T坐标系

径向-横向坐标系指将X-Y坐标系采集得到的地震资料中的每个接收点的X、Y分量旋转一个角度后得到的分量。径向分量(Radial，R)：以震源为起点，并指向接收点方向为R轴正方向，接收P-SV波，随具体的地震道变化。横向分量(Transverse，T)：为R分量旋转90度的方向为T轴正方向，接收P-SH波。

(3)裂缝方位坐标系/自然坐标系/S1-S2坐标系

指转换横波分裂后，快横波偏振方向和慢横波偏振方向。快横波(S1)方向：沿裂缝走向，以与炮检方位角的夹角小于90度的裂缝走向方位角为正向，即S1走向与裂缝走向一致。慢横波(S2)方向：裂缝走向的法线方向，即垂直于裂缝发育方向。

如图1所示，在本发明的一个示例性实施例中基于等级相关分析的多分量转换波裂缝预测方法可以包括：

步骤S01，对转换波三维三分量地震数据的水平分量X和Y进行旋转，得到旋转后的转换波径向分量地震数据和横向分量地震数据，并对转换波进行各向异性处理操作。

在本示例性实施例中，首先在野外进行三维三分量地震数据的采集。在三维三分量地震勘探中，地震采集方法是采集到高质量多分量原始资料的技术保障。由于转换波下行波和上行波具有不对称的射线路径，且转换波的转换点与纵横波速度比和深度有关，因此转换波三维三分量观测系统设计有其自身的特点。

以上，在三维三分量地震勘探中可以使用三分量接收器，三分量接收器分别为水平的X、Y分量接收器，以及垂直的Z分量接收器。由于Inline和Xline方向与横波极化方向不一致，造成X、Y方向上记录的横波既有P-SV波也有P-SH波，没有明确的偏振含义，各个偏振方向上的波场相互混杂，不利于转换波资料的处理及横波分裂信息的研究和提取。因此，为了获得一致的转换波场，在处理三维三分量数据时，需要将水平方向的两个分量能量重新分配，将转换波坐标旋转，实现由野外采集坐标系(X-Y坐标系)向径向-横向坐标系(R-T坐标系)旋转。通过坐标旋转处理后，转换波主要能量分布在径向分量上，横向分量主要代表各向异性的影响。X-Y坐标系向R-T坐标系旋转的示意图如图2所示，在图2中设炮点和检波点连线方向与Inline方向夹角为θ。当径向分量R与X分量方向的夹角0°＜θ＜90°时，坐标旋转式可以为：

其中，T为横向分量。

当径向分量R与X分量方向的夹角90°＜θ＜180°时，坐标旋转式可以为：

R＝-Xcos(180°-θ)+Ysin(180°-θ)，

T＝-Xsin(180°-θ)-Ycos(180°-θ)，

表达成矩阵形式为：

对获取的原始多波数据进行处理，转换波处理是多分量转换波地震勘探能否取得成功的关键。对转换波数据处理，为了得到高信噪比、高保真、无畸变的动校道集，这里，应当注意避免使用道均衡或类似于道均衡之类的均衡手段。多分量转换波数据处理可以包括以下几个关键技术：1、转换波静校正，包括转换波表层静校正和剩余静校正；2、去噪处理；3、地表一致性处理；4、速度分析；5、叠前时间偏移处理等。在数据处理过程中，由于转换横波在近零偏移距时能量较弱，而远偏移距振幅受入射角影响会发生畸变，因此要对偏移后的道集数据进行偏移距选择，以达到较高信噪比的叠加数据。

需要说明的是，转换波静校正、去噪处理、地表一致性处理、速度分析、剩余静校正、叠前时间偏移处理等均是在转换波坐标旋转后进行处理。为了获得高信噪比、高保真、无畸变的动校道集，对转换波数据处理的方法不限于此，其它常规对转换波进行各向异性处理的操作均可。这里的转换波静校正、去噪处理、地表一致性处理、叠前时间偏移处理等均属于常规操作。

步骤S02，选取一个包含目标深度范围的计算时窗，将待处理采样点径向分量地震数据和横向分量地震数据按照计算时窗范围进行截取。

在本示例性实施例中，计算时窗的选取需包含目标深度范围，即包含有裂缝性各向异性地层，可以通过层位控制或者时间段控制目标深度范围。从待处理采样点开始，当计算时窗确定后，将经步骤S01得到的所有采样点的径向分量地震数据和横向分量地震数据分别按照确定的计算时窗范围进行截取，截取的数据分别存入不同的数组中以便进行后续的计算。

以上，待处理采样点的径向分量地震数据和横向分量地震数据是经过处理后的多分量转换波数据。在同一工区中，所有采样点的计算时窗范围选取的层位范围是一致的，例如，通过层位控制，所有采样点的时窗范围均取包含层位1与层位2，但每个采样点对应的范围值在层位1和层位2中是不一样的，比如采样点C的范围可以取1200ms(层位1)～1300ms(层位2)，采样点D的范围可以取1300ms(层位1)～1400ms(层位2)。

步骤S03，对待处理采样点进行等级相关分析，得到等级相关系数最大值。

一方面，由于传统的分析方法多建立在数据变量为定量且服从正态分布的前提下，如果上述前提条件不成立，则会造成分析结果不可信或者错误。另一方面，传统的互相关法计算的是两变量的线性相关，用线性相关描述具有一定的局限性，并且传统的互相关法在抗噪性表现一般。因此，在本发明方法的地震数据分析过程中引入了等级相关分析进行计算，等级相关分析是一种非参数统计相关的分析方法，是一种计算的波形间的相似性是非线性的分析方法。

在本示例性实施例中，所述等级相关分析可以为斯皮尔曼等级相关分析。斯皮尔等级相关也可以称为斯皮尔等级相关系数，是测定变量之间等级相关程度的一种非参数统计相关分析方法。

斯皮尔曼等级相关的基本思路可以是：假设x，y是抽自两个不同总体X，Y的样本，其观察值分别为x₁,x₂,x₃,...,x_n和y₁,y₂,y₃,...,y_n，将他们配对形成(x₁,y₁)，(x₂,y₂)，……(x_n,y_n)。将x_i和y_i各自排序(同时为升序或者降序)，分别计算出它们的秩，计做P_i和Q_i，得到n对秩(P₁,Q₁),(P₂,Q₂),……(P_n,Q_n)。当X和Y完全相关时，∑(P_i-Q_i)＝0，可以记做∑D_i＝0，其中D_i就可以度量X和Y的相关性，若D_i越大，X和Y之间相关越不完全。然而D_i的值可以为正值，也可以为负值，直接用∑D_i来测量相关性会缩小P_i和Q_i之间的差值，因此可以使用∑D_i²来表示P_i和Q_i的差值大小。然而∑D_i²既受到P_i和Q_i不一致程度的影响，同时也受到X和Y样本个数n的影响。鉴于此，为了更加准确地进行相关性分析，采用∑D_i²的最大值去除∑D_i²，得到一个相对的测量指标，叫做斯皮尔曼等级相关，用r_s表示，其中，

其中和分别是P_i和Q_i的平均值。

当每个数据的秩次均不相同时，上式可以写为：

其中，D_i＝P_i-Q_i，即为两个观察值秩之间的差距。

在上述式中，r_s的取值范围在-1到1之间。当r_s大于0时，两者为正相关；当r_s小于0时，两者为负相关。r_s值越接近于1，样本间的相关性越大，r_s越接近于0，样本间的相关性越小。

例如，在本示例性实施例中，假设采样点A(A点的线号为X，道号为Y)的多分量地震数据的采样率为S，计算时窗范围采用层位控制方法，以已知层位hor1为时窗顶，hor2为时窗底。设置该采样点处的旋转角度循环步长为Δθ，时移循环步长为Δτ，阈值范围为θ_k≤θ≤π，τ_k≤τ≤τ_T；其中阈值范围作为循环的终止条件，θ为旋转角度，τ为时移，对采样点A进行等级相关分析，得到等级相关系数最大值的步骤可包括：

(A)根据选取的计算时窗，将截取的待处理采样点A的径向分量和横向分量地震数据分别记录到数组R和数组T中；

(B)对数组R和数组T进行旋转角度θ_j和时移τ_j的变换，得到变换后的分量LR(θ_j,t)和LT(θ_j,t+τ_j)，其中，旋转角度变换式为：

LT(θ_j,t+τ_j)＝LT(θ_j,t)*τ_j

其中，θ_j＝θ_k+Δθ，τ_j＝τ_k+Δτ，θ_k为起始旋转角度，τ_k为起始时移，t为时间。

一般而言，起始角度θ_k可以是0°，起始时移τ_k可以是0ms。

(C)针对不同的旋转角度和时移，分别求取LR(θ_j,t)分量和LT(θ_j,t+τ_j)分量的秩，记为ZR_i和ZT_i；

这里，首先对LR(θ_j,t)分量和LT(θ_j,t+τ_j)分量进行编秩操作，然后求两者的等级相关系数。不同的旋转角度和时移对应着不同的分量，因此，当分量不同时，求取的秩不相等，可以得到不同的ZR_i和ZT_i。

旋转角度和时移必须在设定的阈值范围内进行循环取值，当旋转角度或时移超出阈值范围后，将停止等级相关分析计算。

(D)利用式1计算等级相关系数，得到最大值r_max＝max{r(θ_j,τ_j)}，其中，式1为

其中，D_i＝ZR_i-ZT_i，n为分量内数据的个数。

在本发明方法的等级相关计算过程中，r(θ_j,τ_j)值越大说明旋转后的径向分量与横向分量数据越相关，当求得r(θ_j,τ_j)的最大值r_max时，此时对应的旋转角度可以初步判定是裂缝的发育方位。

步骤S04，判断待处理采样点处是否具有裂缝发育，确定最终裂缝发育方位。

为了剔除裂缝不发育区域，确保对裂缝发育进行统计时更为准确，在进行裂缝发育图绘制时使裂缝带的分布和发育情况更加突出，需要对是否具有裂缝发育作进一步精细的判断，以确定最终的裂缝发育方位。

在本示例性实施例中，对步骤S03等级相关分析结果中的最大值r_max做进一步的判断，以确定待处理采样点处是否有裂缝发育。判断方法可以包括：

步骤S04.1，将待处理采样点的最大值r_max与判别值T₁进行比较，若r_max＞T₁，则进入下一步；

步骤S04.2，计算步骤S03中的待处理采样点所有r(θ_j,τ_j)值的标准差D_r，若D_r＞T₂，则判定该采样点处有裂缝发育，并记录最大值r_max下的旋转角度θ_max和时移τ_max。

其中，判别值T₁，T₂可以是经验值，也可以是给定值。

以上，当判定该采样点处有裂缝发育时，此时的最大值r_max对应的旋转角度即为裂缝发育的方位，时移τ_max即为分裂后快慢波的时差。当判定采样点处的裂缝不发育，则将不保存该采样点的数据，不绘制裂缝发育图。所述判别值T₁，T₂为本领域的常规数值。

因此，在本发明方法的裂缝预测过程中，要判别采样点处是否有裂缝发育，不但需要计算得到最大的等级相关系数，同时还需要满足以上的判别条件。

步骤S05，对判定有裂缝发育的待处理采样点进行快横波和慢横波分离。

在本示例性实施例中，径向分量地震数据和横向分量地震数据均包含有快、慢横波信息，因此，需要将快、慢横波从径向分量地震数据和横向分量地震数据中分离出来。

当确认采样点处有裂缝发育时，可以采用三角变换方法或Alford旋转方法。三角变换方法和Alford旋转方法均是通过坐标旋转的方式对裂缝处的横向分量地震数据以及径向分量地震数据进行分离得到快横波和慢横波。

以上，三角变换方法即是采用式得到快横波和慢横波，其中，S₁代表快横波，S₂代表慢横波。

Alford旋转方法的基本思路可以包括：设裂缝走向方位角β与径向分量正向方位角α的夹角为θ＝β-α，

定义正交旋转矩阵R为：

构造接受信号矩阵：

其中，R₁、T₁和R₂、T₂分别为两个不同方位角的径向量和横向量，两个方位角之间满足正交关系。

根据Alford旋转理论：V＝RSR^T

其中，S₁代表快横波，S₂代表慢横波。

对式V＝RSR^T进行正交旋转后，得到S＝R^TVR。

步骤S06，根据快横波和慢横波时差计算方位各向异性系数。

具体来讲，裂缝越发育，转换波分裂时快慢波时差越大，因此，快慢波时差可以作为裂缝发育程度的一种重要指标，用来计算方位各向异性系数。

当判定采样点处有裂缝发育后，即可进行方位各向异性系数的计算。假设分离后的快波、慢波剖面反射波的存在时间分别为T_S1和T_S2，则快波、慢波时差ΔT＝T_S1-T_S2，由此可以确定各向异性系数为：

K_C＝ΔT/T_S1

这里，快波、慢波时差ΔT即为待采样点处的最大值τ_max。

如图3所示，图3表示横波分裂与裂缝发育方位关系图。横波分裂后快横波沿裂缝发育方位，慢横波与裂缝发育方位垂直，根据波的叠加原理，径向分量和横向分量数据中均包含有快慢波分量。

步骤S07，重复步骤S02至步骤S06，循环至下一个待处理采样点至地震数据所有采样点计算完毕。

这里，还需要对待处理采样点的线号和道号进行判断，当满足start_inline＜X＜end_inline，start_crossline＜Y＜end_crossline时，则可以循环到下一个待处理采样点至地震数据计算完毕，其中，X为待处理采样点的线号，Y为待处理采样点的道号，start_inline是工区起始线号，end_inline是工区终止线号，start_crossline是工区起始道号，end_crossline是工区终止道号。

需要说明的是，在循环处理采样点时，采样点循环的顺序可以依据线号和道号进行循环。例如，在循环过程中，先保持某线号X不变，在道号范围start_crossline＜Y＜end_crossline内依照道号依次循环，道号循环完成后，线号变为X+1，再依次循环所有道号。

步骤S08，绘制裂缝发育预测图。

在本示例性实施例中，可以采用多种模式进行裂缝发育的综合绘制以及对井周围不同范围内的裂缝发育情况进行统计分析。例如，综合裂缝发育的方位以及各向异性系数绘制裂缝发育预测图，每一采样点的裂缝发育用短线表示，短线的指示方向代表该采样点的裂缝发育方位。裂缝发育的强度可以选择以下三种方式绘制：①用直线的长短来表示；②将直线赋予一种颜色，透明度高代表裂缝不发育，透明度越低裂缝越发育；③采用不同色彩代表不同的裂缝发育程度，即用长短、透明度、颜色三种不同方式表示。此外，可以对裂缝发育带的结果进行验证，在进行验证时，需要与测井技术中的裂缝发育情况进行对比，因此可以对井周围一定范围内的裂缝发育情况进行统计，制作玫瑰图，方便与测井数据进行对比分析。

综上所述，本发明的一方面将等级相关理论引入快慢波分离计算中，用非线性相关代替了传统的线性相关，规避了传统线性相关方法对地震数据分布的要求，提高了计算精度，更好地描述了快慢波间的相关性；另一方面采用了分析-判别的流程进一步提高了裂缝预测的准确性，通过设置判别值对裂缝发育计算结果做进一步判别，与传统方法相比，提高了裂缝预测的准确度，突出了裂缝发育较强的地方，能够有效地预测断裂带发育；再一方面采用了多种显示方式综合显示了裂缝防御带的方位和强度，更加直观的显示了裂缝发育情况。

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