一种基于最小熵旋转法的裂缝预测方法与流程

本发明涉及地震勘探技术领域，尤其涉及一种基于最小熵旋转法的裂缝预测方法。

背景技术：

当地震横波在含有裂缝的各向异性介质中传播时，会发生横波分裂现象，形成平行于裂缝走向偏振的快横波和垂直于裂缝走向偏振的慢横波。由于快、慢横波在速度上有差异，因此可以利用它们之间的差异对裂缝属性进行识别。快、慢横波的极化方向反应了裂缝的走向，快、慢横波的走时时差反应了裂缝密度的大小，时差越大，裂缝密度越大。因此研究快慢横波分裂现象成为了研究裂缝方向及其发育程度的最直接有效的方法之一。

目前，现有技术中通过横波分裂法确定裂缝参数的方法有很多种。从对数据的要求来看，可以分为多源法和单源法：

1)多源法以Alford旋转为代表，适用于横波震源激发的四分量VSP数据、九分量数据等，其中VSP是指垂直地震剖面；

2)单源法(互相关法；参数反演法，最小熵旋转法，正交基旋转法等)，适用于P波震源激发，三分量接收的VSP数据。由于横波震源激发成本较高，现实中我们接触的VSP数据多为P波震源激发，所以以P波震源激发的转换横波分裂法的研究更具有实用性。

目前，以P波震源激发的最小熵旋转法比较流行。然而，该方法只能计算得到裂缝的方位角，不能直接获得快慢横波的时差，需要再通过其他方法对分裂后的快慢波进一步分析计算才能得到快慢横波分裂时差，而且同时计算得到的裂缝方位角具有多解性。

技术实现要素：

针对上述问题，本发明提出了一种新改进的基于最小熵旋转法的裂缝预测方 法。该方法包括以下步骤：

S10、对垂直地震剖面数据的水平分量X和Y进行旋转，得到旋转后的径向分量R和横向分量T；

S20、对旋转后的分量进行波场分离，得到分离后的下行分量，并于下行分量上拾取一组发生横波分裂的转换波场。

S30、选取一个目标深度处的检波点，对该检波点接收的每一个震源数据执行以下步骤，以得到该目标深度的裂缝方位角和密度参数：

S30.1、在下行分量上取一个包括了步骤S20所拾取的转换波的时窗，对该时窗内的分量数据进行快慢波分离，基于改进的最小熵旋转计算式扫描计算方位-时差谱，

S30.2、从方位-时差谱中拾取裂缝相对测线的方位和快慢横波时差，并计算裂缝的自然方位角和密度；

S40、对于其他目标深度处的检波点，重复步骤S30，直到得到所有目标深度的裂缝方位和密度参数。

根据本发明的实施例，上述改进的最小熵旋转计算式为

式中，LR(α，t_n)是扫描能量，t_n是快慢横波时差，S1(α,k)是快横波记录，S2(α,k)是慢横波记录，w_n是时窗大小，k为时窗内的采样点，α是方位角。

根据本发明的实施例，上述步骤S30.1包括以下小步骤：

①在下行分量上选取一个相同的时窗，时窗大小要求能包括步骤S20所拾取的转换波；

②扫描一个方位角α，通过式(1)对时窗内的下行分量数据进行快慢横波分离，计算快横波记录和慢横波记录；

式中，S1(α,t)是快横波记录，S2(α,t)是慢横波记录；t为时窗内样点采样时间；R(t)是径向分量数据，T(t)是横向分量数据；

③扫描一个时差t_n，通过所述改进的最小熵旋转计算式计算并保存扫描能量LR(α，t_n)

式中，LR(α，t_n)是扫描能量，t_n是快慢横波时差，S1(α,k)是快横波记录，S2(α,k)是慢横波记录，w_n是时窗大小，k为时窗内的采样点；

④令t_n＝t_n+Δt，重复步骤③，直至t_n＝w_n结束；其中Δt为时差扫描步长；

⑤令α＝α+Δα，重复步骤②～④，直至α＝π/2结束；其中Δα为角度扫描步长；

⑥根据计算的每个LR(α，t_n)制作成方位-时差能量谱图。

根据本发明的实施例，上述方位角α的取值范围为-π/2≤α＜π/2。

根据本发明的实施例，上述时差t_n的取值范围为0≤t_n＜w_n。

根据本发明的实施例，上述步骤S30.2包括以下小步骤：

①从方位-时差能量谱图上拾取能量最大的值，能量最大值所对应的角度就是裂缝相对测线的方位角α，时差t_n就是快慢横波时差；

②计算裂缝的自然方位角β＝θ+α，其中β∈[0，π)、θ是炮点和检波点连线的方位角，θ∈[0，2π)；

③按照下面式(3)计算裂缝密度e：

e≈γ/1.1 (3)

其中，t_s1和t_s2为一裂缝底层界面拾取的快、慢横波的旅行时；t_n为快慢横波时差。

根据本发明的实施例，可以以正东方向为0°。

根据本发明的实施例，可以以逆时针旋转方向为方位角增加的方向。

根据本发明的实施例，上述步骤S30进一步包括步骤S30.3，利用图表统计该目标深度的检波点的所有裂缝方位和密度参数的分布情况。

根据本发明的实施例，可以利用玫瑰图统计该目标深度的检波点的所有裂缝方位分布；利用条形图统计该目标深度的检波点的所有密度参数分布。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：该方法基于横波分裂裂缝预测原理，对现有的最小熵旋转法横波分裂裂缝预测方法进行改进，能够同时得到裂缝方位角和快慢横波的分裂时差，进而得到衡量裂缝发育程度的裂缝密度，同时还解决了原有的裂缝方位角多解性问题。该方法在VSP裂缝的预测中更具有便利性和实用性。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是本发明实施例中采用的裂缝预测方法的工作流程图；

图2a是裂缝的自然方位与测线方位之间的关系的示意图；

图2b是裂缝的自然方位与相对方位之间的关系的示意图；

图3是本发明实施例中分析的一个裂缝介质模型的示意图；

图4是本发明实施例中方位各向异性正演模拟采集方案的示意图；

图5a～5c是本发明实施例中模拟的测线方位角为0°时的三个分量记录X、Y、Z的示意图；

图6a～6b是本发明实施例中分离的下行R、T分量波场；

图7a～7b是本发明实施例中600米目标深度处转换波所在时窗的R、T分量共检波点道集记录；

图8是本发明实施例中600米目标深度处原始R、T分量的偏振曲线图；

图9是本发明实施例中600米目标深度处通过改进的最小熵旋转计算式计算得到的角度-时差谱；

图10是本发明实施例中600米目标深度处横波分裂后的S1、S2偏振曲线图；

图11a～11b是本发明实施例中600米目标深度处横波分裂后的纯快慢横波记录；

图12是本发明实施例中600米目标深度处于各个方位记录计算的裂缝方位、快慢横波延迟时和密度参数的列表；

图13是本发明实施例中获得600米目标深度处的裂缝方位的统计玫瑰图；

图14是本发明实施例中获得600米目标深度处的裂缝密度的统计条形图；

图15是本发明实施例中获得1900米目标深度处的裂缝方位的统计玫瑰图；

图16是本发明实施例中获得1900米目标深度处的裂缝密度的统计条形图。

图17是本发明实施例中利用改进的最小熵旋转法计算式获得的计算结果。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，以下结合附图对本发明作进一步地详细说明。

图1是本发明一实施例中采用的裂缝预测方法的工作流程图。需要说明的是，该图包含了裂缝预测方法的主要步骤，这些步骤可以根据实施时的具体要求进行增减。

S110，对VSP数据的水平分量X和Y进行旋转，得到旋转后的径向分量R分量和横向分量T(以下也可简称R分量和T分量)。

S120，对R分量和T分量进行波场分离，得到分离后的下行R分量和T分量，并在下行R分量或T分量上拾取一组发生横波分裂的转换波场。

S130，对于一个目标深度处的检波点接收的某一炮数据，在下行R分量和T分量上取一个包括了步骤S120所拾取的转换波的时窗，对该时窗内的数据利用改进后的最小熵旋转法扫描计算角度-时差谱。

S140，从角度-时差谱中拾取裂缝相对测线的方位和快慢波时差，并计算裂缝自然方位角和密度。

S150，对步骤S130的目标深度处检波点接收的其他炮的数据重复执行步骤S130和S140。

S160，根据经过步骤S150后获得的某一目标深度处所有的裂缝方位角绘制玫瑰图，对裂缝密度绘制条形图，以完成方位角和密度参数的统计。

S170，对于其他目标深度处的检波点，重复步骤S130～160，直到得到所有目标深度的裂缝方位和密度参数。

下面对各步骤进行详细说明。其中：

步骤110和120属于地震技术领域中对垂直地震剖面数据的常规处理，此处不再叙述具体做法。

步骤S130进一步包括：

①在下行分量上选取一个相同的时窗，时窗大小w_n要完全包括了步骤S200所拾取的转换波。

②扫描一个方位角α，-π/2≤α＜π/2，通常从-π/2开始扫描。然后，通过下面的式(1)对时窗内的下行分量数据进行快慢横波分离，计算快横波记录S1(α,t)和慢横波记录S2(α,t)。

式中的t为时窗内样点采样时间；R(t)是径向分量数据，T(t)是横向分量数据。

③扫描一个时差t_n，0≤t_n＜w_n，通常从0开始扫描。然后，通过下面改进的最小熵旋转计算式(2)计算扫描能量LR(α，t_n)，并保存结果。

式中，LR(α，t_n)是扫描能量，t_n是快慢横波时差，S1(α,k)是快横波记录，S2(α,k)是慢横波记录，w_n是时窗大小，k为时窗内的采样点。

④令t_n＝t_n+Δt(Δt为时差扫描步长)，重复步骤③，直至t_n＝w_n结束。

⑤令α＝α+Δα(Δα为角度扫描步长)，重复步骤②～④，直至α＝π/2结束。

⑥将步骤③计算的每个LR(α，t_n)按照能量强弱赋予不同颜色深度，制作成方位-时差能量谱图(如图10所示)。当然，方位-时差能量谱图也可以不局限于这种表现方式。

步骤S140进一步包括：

①从方位-时差能量谱图上拾取能量最大的值，能量最大值所对应的角度就是裂缝相对测线的方位角α，时差t_n就是快慢横波时差。

②计算裂缝的自然方位角β＝θ+α，其中β∈[0，π)、θ∈[0，2π)是炮点和检波点连线(R分量)的方位角。为方便计算分析，在本实施例中以正东方向为 0°，以逆时针旋转方向为方位角增加的方向(如图3所示)。

③按照下面式(3)计算裂缝密度：

e≈γ/1.1 (3)

其中，t_s1和t_s2为一裂缝底层界面拾取的快、慢横波的旅行时；t_n为快慢横波时差。

第一个实施例

图2a～2b显示了裂缝的自然方位、相对方位与测线方位之间的关系的示意图。图2a中显示了通过水平分量旋转得到的径向分量R、横向分量T与测线方位之间的关系；图2b中显示了裂缝自然方位与相对方位之间的关系。在图中，N、E分别表示北向和东向。

下面以图3所示的一个裂缝介质模型为例来说明本发明的裂缝预测方法。该模型共有四层介质，其中第二层、第四层为HTI裂缝各向异性层，模型及介质参数如图3所示。图中，Vp是纵波速度，Vs是横波速度，ρ是密度，ε、γ和δ是Thomsen参数。

图4是方位各向异性正演模拟采集方案的示意图。从该图可知，以井口为中心，每隔15°为一个方位，在共24个方位上部署炮点，炮点距井口的偏移距为500米；接收点分布于井中深度100～2000米处，每隔10米放置一个。

图5a～5c是模拟的测线方位角为0°(正东方向，逆时针)的三个分量X、Y、Z记录。其中，在X、Z分量上各个层都观察到了P波；在Y分量上，在各向同性层(1、3层)接收不到P波，而在各向异性层(2、4)中接收到了P波。总之，在X、Y、Z分量都观察到了由震源激发的P波在地层分界面上转换的S波(以下简称转换波)，并且在各向异性层(2、4层)横波发生了分裂现象，故而，对于上行波来说，三个分量上都有反射P波和转换波投影，其中，Y分量上以上行转换波为主，Z分量上以上行P波为主。

图6a～6b是分离的下行R、T分量波场。在该波场中包括了下行P波、第二层HTI转换波、第三层各向同性层转换波和第四层HTI层转换波。其中，虚线标注的是第二层HTI层转换波，这组转换波是准备分析的横波。

为了验证本发明所提出的裂缝预测方法的正确性，这里分别分析了第二层HTI介质600米目标深度检波点的裂缝参数(方位和密度参数)和第四层HTI介 质1900米深度点的裂缝参数。

图7a～7b是600米目标深度处转换波所在时窗的R、T分量共检波点道集记录。图7a中的R分量和图7b中的T分量道集均呈“正弦”变换，道集中的记录从0°至245°。其中，走时在平行裂缝时(50°和230°)最小；R分量能量在平行裂缝时最大，在垂直裂缝时(140°和320°)最小；T分量在平行裂缝和垂直裂缝时极性都发生了一次变化。

图8是600米目标深度处原始R、T分量的偏振曲线图(横轴是R分量，竖轴是T分量)。从图中可以看到，每个偏振图有两组近乎线性偏振且近似垂直的两个偏振方向。其中，能量较大的是快横波偏振方向，能量较小的为慢横波偏振方向。

图9是600米目标深度处通过改进的最小熵旋转计算式计算得到的角度-时差谱(横轴是角度，竖轴是时差)。在该角度-时差谱上能量最大值对应于裂缝与测线夹角方向。

图10是600米目标深度处横波分裂后的S1、S2偏振曲线图(横轴是S1波，竖轴是S2波)。经过横波分裂后，快横波S1被分解到X轴方向，慢横波S2被分解到Y轴方向，快慢横波的偏振方向近似呈正交。

图11a～11b是600米目标深度处横波分裂后的纯快慢横波记录。可以看到，分裂后快慢波走时形态基本一致，波形也相似，同时快横波之间有明显的到时差。

图12是600米目标深度处于各个方位记录计算的裂缝方位、快慢横波延迟时和密度参数的列表。通过对这个表中的数据进行统计，最终获得图13所示的裂缝方位的统计玫瑰图和图14所示的裂缝密度的统计条形图。其中，从图13的方位玫瑰图中可以看出，裂缝方位大致都指向50°，经过统计裂缝方位的平均值为50.08°。这说明预测的裂缝方位是正确的，而且精度较高。从图14的密度条形图中可以得出，裂缝密度的平均值为0.16，基本接近于理论计算的裂缝密度0.17。这说明预测的裂缝密度是正确的，而且精度较高。

采用同样的方法流程，分析了第四层HTI层1900米深度处的裂缝参数。附图15是预测的裂缝方位角玫瑰图，从以上图上看到，裂缝方位角指向110°左右，经过统计裂缝方位平均值为109.96°。附图16是计算的裂缝密度统计图,裂缝密度均值为0.11，和真实裂缝密度0.09基本接近。

上述实施例有力地验证了本发明提出的裂缝预测方法是有效的。该方法弥补了现有技术中最小熵旋转法不能直接得到快慢横波分裂时差和方位具有多解性的不足。改进的最小熵旋转计算式的右端的前半部分是未改进的最小熵旋转公式的简化形式：

其只与角度有关，与时差无关；而改进的最小熵旋转计算式既与角度有关，又与时差有关。如图17所示，利用改进的最小熵旋转计算式得到的结果是一条正弦曲线，有两组极大值的解(极大值对应着最小熵)，并且它们之间相差90°。虽然还需要通过其他信息进一步判断其真解(从两组极大值的解选择一个)，但是已经能够改进传统方法中方位角具有多解性的不足，对实际工程勘探有很好的指导意义。

以上所述，仅为本发明的具体实施案例，本发明的保护范围并不局限于此，任何熟悉本技术的技术人员在本发明所述的技术规范内，对本发明的修改或替换，都应在本发明的保护范围之内。