多分量转换波裂缝预测方法与流程

本发明属于石油勘探领域，涉及一种针对裂缝性油气储层的裂缝系统检测方法，更具体地讲，涉及一种多分量转换波裂缝预测方法。

背景技术：

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

随着多波研究的不断深入，多波勘探技术越来越多地应用到油气勘探开发中。横波分裂作为天然地震、石油勘探领域的一项具有强大潜力的技术被地质学家和地球物理学家广泛关注。横波分裂，又称为横波双折射，即当横波穿过方位各向异性介质时会分裂成沿平行于裂隙方向传播的快横波和沿垂直裂隙方向传播的慢横波，两者传播速度不同，因此在各向异性介质中形成快慢横波时差，振幅衰减等，这些特性都与各向异性性质有关。可根据快横波、慢横波的时差、波形、振幅衰减等反过来研究裂缝系统的方位和密度。利用横波分裂预测储层裂缝成为一种直接可靠的方法。但是横波震源比较昂贵，利用纯横波震源来进行横波分裂勘探还难以实现大规模应用。转换波勘探克服了纯横波勘探激发难、成本高、静校正量大等缺陷，同时具有信息丰富、兼有纵波和横波的优势、资料信噪比较高等优点，使得多分量转换波勘探成为油气储层探测的有力工具。

在转换波勘探裂缝预测中，准确地计算出裂缝发育方位和时差对于准确刻画裂缝系统在地下的分布形态具有重要意义。目前，常用的方法有互相关法和能量比法。传统方法中计算的相关性或者协方差矩阵时，均采用线性相关的方式，而波形间的相似性是非线性的，且参数计算的数据样本容量有限，不满足正态分布的特征，从数学意义上讲，样本数据本身就不满足要求。

技术实现要素：

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

为了实现上述目的，本发明的提供了一种多分量转换波裂缝预测方法，所述预测方法可以包括以下步骤：选取一个包含目标深度范围的分析窗口，将待处理采样点的径向分量地震数据和横向分量地震数据按照分析窗口进行截取；对截取的径向分量地震数据和横向分量地震数据进行角度旋转，并对角度旋转后的横向分量地震数据进行时延转换；构建目标函数，计算目标函数最小值；记录目标函数最小值对应的旋转角度和时延；对待处理采样点进行快横波和慢横波分离，根据快横波和慢横波的旅行时差计算裂缝发育强度；重复所述选取分析窗口的步骤至所述计算裂缝发育强度的步骤，循环至下一个待处理采样点直到地震数据计算完毕，其中，所述构建目标函数，计算目标函数最小值的步骤包括：

记截取的径向分量地震数据进行角度旋转后的分量为p(θi+n,t)，角度旋转后的横向分量地震数据进行时延转换的分量为q(θi+n,t+ti+n)；分别对p(θi+n,t)、q(θi+n,t+ti+n)进行编秩操作，记为rgp和rgq，其中，θi+n表示旋转角度，t表示时间，ti+n表示时延，n表示整数；构建矩阵其中，m表示p或q，n表示p或q,ρrgm,rgn表示皮尔森相关系数，cov(rgm,rgn)表示求协方差函数，σrgp表示rgp的标准差，σrgq表示rgq的标准差；计算矩阵u(θi+n,t+ti+n)的特征值eig(u(θi+n,t+ti+n))；其中，θi+n＝θi+nδθ，ti+n＝ti+nδt，θmin≤θi+n≤θmax且tmin≤ti+n≤tmax，θi为起始旋转角度，ti为起始时延，δθ为旋转角度循环步长，δt为时延循环步长，[θmin,θmax]为旋转角度阈值范围，[tmin,tmax]为时延阈值范围；构建目标函数f(θ,t)＝min[eig(u(θi+n,t+ti+n))]，n从0开始取值，n每次增加1，得到目标函数最小值，其中，min[eig(u(θi+n,t+ti+n))]表示求特征值eig(u(θi+n,t+ti+n))的最小值。

与现有技术相比，本发明将非参数统计理论引入多波裂缝预测过程中，采用等级相关的方式替代传统算法中的线性相关方法，等级相关方式对于波形相似性的测度更加敏感，同时规避了线性相关对样本数据的要求，提高了计算的准确度，为裂缝性油气储层的精细描述提供依据。

附图说明

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

图1示出了根据本发明示例性实施例的多分量转换波裂缝预测方法流程示意图。

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

具体实施方式

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

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

本发明的创新之处在于，本发明的方法将非参数统计的方法引入裂缝预测方法中，在构建目标函数时对转换后的数据进行编秩操作，然后构建基于秩的协方差矩阵，将线性相关变为非线性相关，提高了计算的准确度，为地下裂缝识别和预测提供更可靠依据。

图1示出了根据本发明示例性实施例的多分量转换波裂缝预测方法流程示意图。图2示出了根据本发明示例性实施例的横波分裂与裂缝发育方位关系示意图。

本发明提供了一种多分量转换波裂缝预测方法，如图1所示，在本发明的一个示例性实施例中，所述方法可以包括：

步骤s01，对转换波三维三分量地震数据的水平分量x和y进行旋转，得到旋转后的转换波径向分量地震数据r(即径向分量r)和横向分量地震数据t(即横向分量t)。

以上，在三维三分量地震勘探中可以使用三分量接收器，三分量接收器分别为水平的x、y分量接收器，以及垂直的z分量接收器。由于inline和crossline方向与横波极化方向不一致，造成x、y方向上记录的横波既有p-sv波也有p-sh波，没有明确的偏振含义，各个偏振方向上的波场相互混杂，不利于转换波资料的处理及横波分裂信息的研究和提取。因此，为了获得一致的转换波场，在处理三维三分量数据时，需要将水平方向的两个分量能量重新分配，将转换波坐标旋转，实现由野外采集坐标系(x-y坐标系)向径向-横向坐标系(r-t坐标系)旋转。通过坐标旋转处理后，转换波主要能量分布在径向分量上，横向分量主要代表各向异性的影响。坐标旋转式可以为：

其中，径向分量r与x分量方向的夹角或

可以对旋转后的转换波进行各向异性处理操作，所述各向异性处理操作属于本领域的常规方法。

步骤s02，确定分析窗口范围，根据分析窗口截取待处理采样点的径向分量地震数据r和横向分量地震数据t。

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

以上，在同一工区中，所有采样点的分析窗口范围选取的层位范围是一致的，例如，通过层位控制，所有采样点的时窗范围均取包含层位1与层位2，但每个采样点对应的范围值在层位1和层位2中是不一样的，比如采样点c的范围可以取1200ms(层位1)～1300ms(层位2)，采样点d的范围可以取1300ms(层位1)～1400ms(层位2)。

步骤s03，对径向分量地震数据和横向分量地震数据进行旋转角度和时延转换，得到转换后的分量。

在本示例性实施例中，由于径向分量和横向分量中记录的数据是快横波和慢横波信息在地表接收时的重新组合，即径向分量和横向分量中均为快横波和慢横波的信息的叠加。因此，需要通过旋转的方式将快横波和慢横波进行分离。同时快横波和慢横波均是由同一个转换横波分裂而得，它们的波形是相似的，慢横波的传播速度比快横波传播速度慢，两者之间有时差。因此，在角度旋转的基础上，需要对慢横波进行时延补偿。

具体的，首先分别对待处理采样点处的截取的径向分量r和横向分量t进行角度旋转，旋转式可以为：

得到分量p(θi+n,t)和q(θi+n,t)，然后对旋转后的横向分量数据进行时延，q(θi+n,t+ti+n)为q(θi+n,t)中的每个值向上时移ti+n/s个点，s为采样率，θi+n为旋转角度，ti+n为时延，θi为起始旋转角度，ti为起始时延，n表示整数，θi+n＝θi+nδθ，ti+n＝ti+nδt，θmin≤θi+n≤θmax且tmin≤ti+n≤tmax，δθ为旋转角度循环步长，δt为时延循环步长，[θmin,θmax]为旋转角度阈值范围，[tmin,tmax]为时延阈值范围。

以上，q(θi+n,t+ti+n)为q(θi+n,t)中的每个值向上时移ti+n/s个点指的是，例如，时移ti+n/s取1个点，则q(θi+n,t)中的数据组q1,q2,q3,q4...时移一个点，q(θi+n,t+ti+n)对应的数组变为从q2开始，即为q2,q3,q4...。

从起始角度和起始时延开始，即起始值开始，通过角度旋转式和时延计算式，得到截取的径向分量地震数据角度旋转后的分量为p(θi+n,t)，并对角度旋转后的横向分量地震数据进行时延转换后的分量为q(θi+n,t+ti+n)。

这里，时延计算主要针对旋转后的横向分量，即对角度旋转后的横向分量地震数据p(θi+n,t)进行时延转换后的分量为q(θi+n,t+ti+n)。旋转角度θi+n和时延ti+n均需要在阈值范围内取值，阈值范围可以是设定值或经验值。例如，旋转角度的阈值范围一般可以为[-π,π]或[0,2π]。

例如，设处理后的多分量数据的采样率为s，分析窗口通过层位控制，层位hor1为时窗顶，hor2为时窗底。从采样点a(p,q)开始循环计算，其中，p表示采样点a的线号，q表示采样点a的道号，将径向分量数据和横向分量数据按照上述分析窗口范围进行截取，截取的数据分别存入数组r＝[r1,r2,r3,......,rn]和数组t＝[t1,t2,t3,......,tn]中进行计算。

然后，对数组r和t进行θi+n角度的旋转，例如，可以设旋转角度的阈值范围为0＜θi+n＜π；然后对旋转后的横向分量数据进行时延补偿，时延的阈值范围可以设为[0,t]，得到截取的径向分量地震数据角度旋转后的分量为p(θi+n,t)，角度旋转后的横向分量地震数据进行时延转换后的分量为q(θi+n,t+ti+n)。

步骤s04，构建目标函数，求取目标函数最小值。

在本示例性实施例中，目标函数的构建可以包括：

步骤s04.1，首先，对分量p(θi+n,t)和q(θi+n,t+ti+n)进行编秩操作，即分别求取分量p(θi+n,t)和q(θi+n,t+ti+n)的秩，记为rgp和rgq，其中，θi+n表示旋转角度，t表示时间，ti+n表示时延，n表示整数。

步骤s04.2，构建矩阵其中，

m表示p或q，n表示p或q,即

ρrgm,rgn表示皮尔森相关系数，cov(rgm,rgn)表示求协方差函数，σrgp表示rgp的标准差，σrgq表示rgq的标准差。

步骤s04.3，计算矩阵u(θi+n,t+ti+n)的特征值eig(u(θi+n,t+ti+n))。其中，θi+n＝θi+nδθ，ti+n＝ti+nδt，θmin≤θi+n≤θmax且tmin≤ti+n≤tmax，θi为起始旋转角度，ti为起始时延，δθ为旋转角度循环步长，δt为时延循环步长，[θmin,θmax]为旋转角度阈值范围，[tmin,tmax]为时延阈值范围。

这里，需要判断旋转角度θi+n和时延ti+n是否在阈值范围内，如果满足阈值条件，则可以继续循环取新的旋转角度和时延值，带入步骤s03，得到不同的分量p(θi+n,t)和q(θi+n,t+ti+n)，进而得到不同的特征值。

步骤s04.4，构建目标函数f(θ,t)＝min[eig(u(θi+n,t+ti+n))]，n从0开始取值，n每次增加1，得到目标函数最小值，其中，min[eig(u(θi+n,t+ti+n))]表示求特征值eig(u(θi+n,t+ti+n))的最小值。

步骤s05，记录目标函数f(θ,t)取最小值时对应的旋转角度和时移。

这里，当目标函数取最小值时，此时的旋转角度即为裂缝发育的方位角，时延即为快慢波的时差。转化波分裂时，快慢波时差越大，表明裂缝越发育，因此，快慢波的时差可以作为裂缝发育程度的一种重要指标。

步骤s06，对待处理采样点进行快横波和慢横波分离，根据快横波和慢横波的旅行时差计算裂缝发育强度。

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

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

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

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

定义正交旋转矩阵g为：

构造接受信号矩阵：

其中，r1、t1和r2、t2分别为两个不同方位角的径向量和横向量，两个方位角之间满足正交关系。

根据alford旋转理论：v＝rsr^t，

其中，s1代表快横波，s2代表慢横波，

对式v＝rsr^t进行正交旋转后，得到s＝r^tvr。

所述计算裂缝发育强度的方法包括：裂缝发育强度kc＝δt/ts1，其中，δt＝ts1-ts2，ts1为分离后的快横波剖面反射波存在时间，ts2为分离后的慢横波剖面反射波存在时间。

步骤s07，重复步骤s01至步骤s06，循环至下一个待处理采样点直到地震数据计算完毕。

这里，在处理下一个待采样点时，还需要对上一个处理的采样点的线号和道号进行判断，例如，假设采样点a(p,q)，当满足start_inline＜p＜end_inline，start_crossline＜q＜end_crossline时，则可以循环到下一个待处理采样点，直到地震数据计算完毕，其中，p为待处理采样点的线号，q为待处理采样点的道号，start_inline是工区起始线号，end_inline是工区终止线号，start_crossline是工区起始道号，end_crossline是工区终止道号。

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

步骤s08，绘制裂缝发育综合分布图。

这里，在所有采样点计算完毕后，对每个采样点的数据结果进行保存，形成两个数据，即裂缝发育方位角和裂缝发育强度，然后在综合裂缝发育方位和裂缝发育强度绘制裂缝发育综合分布图。所述绘制裂缝发育综合分布图的方法包括：每一采样点的裂缝发育用短线表示，短线的方位表示裂缝发育的方位角，短线的长短表示裂缝发育强度。

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

综上所述，本发明通过将非参数统计方法引入到多分量转换波裂缝预测中，规避了传统线性相关方法对地震数据分布的要求，在进行协方差矩阵计算时，代替了等级相关计算结果，能够很好地描述两矢量间的相关性，提高了计算准确度，使预测结果更加准确、可靠。

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