CN115903039B

CN115903039B - 一种基于进退法的横波分裂分析方法

Info

Publication number: CN115903039B
Application number: CN202211677314.6A
Authority: CN
Inventors: 周怀来; 杨宇勇; 何正涛; 王元君
Original assignee: Chengdu Univeristy of Technology
Current assignee: Chengdu Univeristy of Technology
Priority date: 2022-12-26
Filing date: 2022-12-26
Publication date: 2023-10-24
Anticipated expiration: 2042-12-26
Also published as: CN115903039A

Abstract

本发明公开了一种基于进退法的横波分裂分析方法，主要步骤：首先、根据快慢横波相关性构造目标函数；然后，通过确定的真实值(或实际数据)得到R分量数据和T分量数据，进一步计算得到S1(t)和S2(t)；根据区间大小，将区间划分为4个子区间，取每个子区间的中心点作为初始点；在每个区间采用进退法确定该区间的结果值，得到4个结果值，比较4个结果值，取最大COV值对应的时差和角度作为输出值。本发明以快慢波的相关性作为目标函数，通过缩小极值区间来确定极值点，并且通过分区的方式获得真实的裂缝参数，为非常规油气储层提供准确的裂缝信息。

Description

一种基于进退法的横波分裂分析方法

技术领域

本发明涉及多分量地震资料处理技术领域，尤其是一种针对裂缝型储层预测的横波分裂分析方法。

背景技术

裂缝作为检测非常规油气含量的重要参数，它不仅可以连接孤立的孔隙，增加储层的有效孔隙度，还可以作为重要的油气运移空间，提高储层的渗透率。20世纪90年代开始，随着多波多分量勘探技术开始登上舞台，横波分裂技术开始应用转换横波，近年来，随着3D3C多波多分量勘探技术的飞速发展，横波分裂技术的发展又到了一个新的高潮。当横波的偏振方向以斜交地方式穿过含裂缝的地层时，会发生横波分裂现象(SWS)。初始横波分裂为一个快横波和一个慢横波，分裂的快横波偏振方向与裂缝走向平行，快慢波时差与裂缝的强度密切相关，因此，横波分裂裂缝检测技术成为研究裂缝方向和其发育程度最直接最可靠的方法之一。当今常规的横波分裂分析方法如旋转相关法、切向能量法、协方差矩阵法都是一种枚举式的网格扫描算法。旋转相关法通过快慢横波相关性能很好的寻找最佳的裂缝参数(快横波偏振方向和快慢横波时差)。但是，在面对庞大的三分量采集数据时，枚举式的计算特点必定会造成计算效率低下。所以，如何快速地提高计算效率，得到准确的裂缝参数信息，对裂缝型油气藏储层预测十分重要。

发明内容

针对现有的横波分裂分析方法存在的面对庞大的三分量采集数据时计算效率低下的问题，本发明提供一种计算效率快的新的基于进退法的横波分裂分析方法，该方法以快慢波的相关性作为目标函数，通过缩小极值区间来确定极值点，快速而准确的找到真实的裂缝参数。从而，得到计算速率最快，精度最高的横波分裂分析方法，为非常规油气储层提供准确的裂缝信息。

本发明提供的基于进退法的横波分裂分析方法，步骤如下：

S1、根据快慢横波相关性构造目标函数如下：

式中，θ和Δt为SWS(横波分裂)参数，S1(t)和S2(t)分别为一个时窗内的快横波和慢横波的时间序列；和/>分别是S1(t)和S2(t)在一个时窗内的平均值。

S2、通过确定的真实值或实际数据得到此时的R分量数据R(t)和T分量数据T(t)，由R(t)和T(t)根据以下公式计算得到S1(t)和S2(t)：

S3、将测试区间划分为4个子区间，取每个子区间的中心点作为初始点。例如总区间为180×60，180为方位角范围，60为时差范围，分为四个子区间后，每个子区域的大小为90×30。

S4、根据给定的初始点确定区间，采用二维进退法确定该区间的结果值。具体方法步骤如下：

(1)对Δt方向零点的搜索；

第一步、确定范围，范围的上下界由初始点(θ_k，Δt_k)和比较点(θ_k，Δt_k+1)确定；Δt_k+1等于Δt_k加上步长h；如果COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)小于零，得到范围[Δt_k，Δt_k+1]，进行下一步；如果COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)大于零，得到新的Δt_k+1，新Δt_k+1等于新的Δt_k加上新的h，新的Δt_k等于上一步的Δt_k+1，新的h是上一步h的两倍；然后，再次判断COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)的值是大于零还是小于零，确实是执行下一步或者继续更新这两个点；

第二步、确定另外两个比较点(θ_k，l_t)和(θ_k，u_t)缩小范围；其中，l_t和u_t的计算公式如下：

l_t＝Δt_k+0.382×(Δt_k+1-Δt_k)，

u_t＝Δt_k+0.618×(Δt_k+1-Δt_k).

通过比较|COV′_Δt(θ_k，l_t)|与|COV′_Δt(θ_k，u_t)|的大小，更新lt或者u_t：

如果|COV′_Δt(θ_k，l_t)|>|COV′_Δt(θ_k，u_t)|，则

Δt_k＝l_t，

l_t+1＝Δt_k+0.382×(Δt_k+1-Δt_k)

如果|COV′_Δt(θ_k，l_t)|<|COV′_Δt(θ_k，u_t)|，则

Δt_k+1＝u_t，

u_t+1＝Δt_k+0.618×(Δt_k+1-Δr_k)；

继续迭代，直到|COV′_Δt(θ_k，l_t)|与|COV′_Δt(θ_k，u_t)|之间的差值小于精度εΔt，εΔt＝10^-6，即得到Δt方向的零点，Δt＝(l_t+u_t)/2；

(2)对θ方向零点的搜索；

第一步、由初始点(θ_k，Δt)得到比较点(θ_k+1，Δt)，θ_k+1等于θ_k加上步长h；然后，计算COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值，如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)<0，则范围为[θ_k，θ_k+1]，执行下一步；如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)>0，继续更新，得到新的θ_k+1，新的θ_k+1等于新的θ_k加上新的h，新的θ_k，等于上一个θ_k+1，新的h是上一个h的两倍；然后，再次判断COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值是大于零还是小于零，确实是执行下一步或者继续更新这两个点；

第二步、确定另外两个比较点(l_θ，Δt)和(u_θ，Δt)缩小范围；其中，l_θ和u_θ的计算公式如下：

l_θ＝θ_k+0.382×(θ_k+1-θ_k)

u_θ＝θ_k+0.618×(θ_k+1-θ_k)

然后，通过比较|COV′_θ(l_θ，Δt)|和|COV′_θ(u_θ，Δt)|的大小，更新l_t或者u_t：

如果|COV′_θ(l_θ，Δt)|>|COV′_θ(u_θ，Δt)|，则

θ_k＝l_θ

l_θ+1＝θ_k+0.382×(θ_k+1-θ_k)

如果|COV′_θ(l_θ，Δt)|<|COV′_θ(u_θ，Δt)|，则

θ_k+1＝u_θ

u_θ+1＝θ_k+0.618×(θ_k+1-θ_k)

继续迭代，直到|COV′_θ(l_θ，Δt)|和|COV′_θ(u_θ，Δt)|的差值小于精度ε_θ，ε_θ＝10^-6，即得到θ方向的零点，θ＝(l_θ+u_θ)/2；

(3)计算点(θ，Δt)的COV′，如果COV′的值小于精度ε₃，ε₃＝10^-8，停止循环，将该点作为结果；反之，则继续重复步骤(1)和(2)交替搜索方向，寻找新的Δt和θ方向的零点，直到结果满足条件。

如果一个点在正、负两个方向都可以搜索Δt或θ方向的零点，则计算两个方向的零点，并选择较低的COV值作为更新点。

S5、按照步骤S4的方法，在其他三个子区间重复进行，每个子区间得到一个结果值。比较4个结果值，取得到最大COV值对应的时差和角度作为输出值。

与现有技术相比，本发明的有益之处在于：

(1)本发明以快慢波的相关性作为目标函数，通过缩小极值区间来确定极值点，快速而准确的找到真实的裂缝参数，从而得到一种计算速率最快，精度最高的横波分裂分析方法。

(2)因为基于最优化算法容易陷入局部收敛，本发明将测试区间划分为四个子区域，并在每个子区域内进行计算，得到四个输出值，比较四个输出值选择最大的COV值对应的点为结果点，以排除局部最优解。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1、横波分裂(SWS)示意图。

图2、合成数据的原始波数据图。

图3、基于黄金分割法的进退法流程图。

图4、进退法搜索过程及模拟寻找裂缝参数真实值；(a)一维进退法求最大值范围(Δt_k,Δt_k+1)，(b)通过缩小范围来确定零点(最优点)，(c)左下角为初始点。

图5、是进退示意图和区域划分示意图。

图6、本发明基于进退法的横波分裂分析方法的系统框架图。

图7、测试范围内的点依次作为真实值时的合成数据误差图。

图8、原始数据R(a)和T(b)分量剖面图。

图9、实际数据计算结果图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

如图1所示，图中，地面上的两个箭头分别表示快、慢横波的偏振方向，θ表示裂缝方位角。当横波到达反射界面时，会产生上行转换横波，假如横波的偏振方向与裂缝面倾斜时，转换横波会分裂为一组正交的快、慢横波。快、慢横波的传播方向几乎垂直于地面，且快横波到达地面的速度快于慢横波。三分量检波器的R和T分量可以在地面上接收含快慢横波的混合信号。R分量方向与检波器-炮点测线方向平行，T分量方向与其正交。

当转换横波S(t)斜交地穿过裂缝面时，用下式计算波场：

其中，θ和Δt为SWS参数，S1(t)和S2(t)分别为一个时窗内的快横波和慢横波的时间序列。

R和T分量可以通过含θ和S(t)的向量转换得到:

对该公式乘以θ的逆矩阵得到下述公式：

因此，从现场接收的R和T分量可以得到S1(t)和S2(t-Δt)。最后，根据快慢横波的相关关系构造目标函数：

其中，和/>分别是S1(t)和S2(t)在一个时窗内的平均值。然后计算网格上各点的COV值，确定合理的θ和Δt(COV＝1)作为SWS参数。

图2是合成数据的原始波数据，a为S(t)、b为S1(t)和S2(t)、c为S1(t)和S2(t-Δt)、d为R(t)和T(t)。

本发明中以COV为目标函数，通过进退法确定搜索路径。进退法通过上界值和下界值的异号性确定含零点的范围，黄金分割法通过缩小范围得到最大值(极值点)。二维进退法需要规律地搜索两个方向自变量的零点，得到最大的COV值。通过在一个方向上搜索零点，然后继续从此点在另一个方向上搜索零点。如果该点的COV′值小于ε₃(一个很小的正数)，则循环结束；否则，继续搜索零点，直到结果满足条件(COV′小于ε₃)。本发明定义的搜索模式为Δt-θ。

本发明采用的基于黄金分割法的进退法流程如图3所示；主要步骤如下：

(1)对Δt方向零点的搜索；

第一步、在寻找Δt的最大值时，确定范围。范围的上下界由初始点(θ_k，Δt_k)和比较点(θ_k，Δt_k+1)确定；Δt_k+1等于Δt_k加上步长h；

如果COV′_Δt(θ_k，Δ_tk)COV′_Δt(θ_k，Δt_k+1)<0，得到范围Δt∈[Δt_k，Δt_k+1]，进行下一步；

如果COV′_Δt(θ_k，Δt_k)COV′_Δt(θ_k，Δt_k+1)>0，得到新的Δt_k+1，新Δt_k+1等于新的Δt_k加上新的h，即Δt_k+1＝Δt_k+h，新的Δt_k等于上一步的Δt_k+1，Δt_k＝Δt_k+1，新的h是上一步h的两倍，h→2h；然后，再次判断COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)的值是大于零还是小于零，确实是执行下一步或者继续更新这两个点。

第二步、确定另外两个比较点(θ_k，l_t)和(θ_k，u_t)缩小范围；其中，l_t等于Δt_k加上Δt_k+1与Δt_k之差的0.382倍。u_t等于Δt_k加上Δt_k+1与Δt_k之差的0.618倍。即l_t和u_t的计算公式如下：

l_t＝Δt_k+0.382×(Δt_k+1-Δt_k)，

u_t＝Δt_k+0.618×(Δt_k+1-Δt_k).

Δt_k＝l_t，

l_t+1＝Δt_k+0.382×(Δt_k+1-Δt_k)

如果|COV′_Δt(θ_k，l_t)|<|COV′_Δt(θ_k，u_t)|，则

Δt_k+1＝u_t，

u_t+1＝Δt_k+0.618×(Δt_k+1-Δt_k)；

继续迭代，直到|COV′_Δt(θ_k，l_t)|与|COV′_Δt(θ_k，u_t)|之间的差值小于精度εΔt(εΔ_t≤10^-6)，即得到Δt方向的零点，Δt＝(l_t+u_t)/2。这样就完成了对Δt方向零点的搜索。

(2)对θ方向零点的搜索；

从点(θ_k，Δt)开始搜索θ方向的最大值，θ方向的最大值更新过程类似于寻找Δt的最大值。

第一步、由初始点(θ_k，Δt）得到比较点(θ_k+1，Δt)，θ_k+1等于θ_k加上步长h(一个新的步长)；然后，计算COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值；

如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)<0，则范围为θ∈[θ_k，θ_k+1]，执行下一步；

如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)>0，继续更新，得到新的θ_k+1，新的θ_k+1等于新的θ_k加上新的h，θ_k+1＝θ_k+h，新的θ_k，等于上一个θ_k+1，θ_k＝θ_k+1；新的h是上一个h的两倍，h→2h；然后，再次判断COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值是大于零还是小于零，确实是执行下一步或者继续更新这两个点；

l_θ＝θ_k+0.382×(θ_k+1-θ_k)

u_θ＝θ_k+0.618×(θ_k+1-θ_k)

θ_k＝l_θ

l_θ+1＝θ_k+0.382×(θ_k+1-θ_k)

如果|COV′_θ(l_θ，Δt)|<|COV′_θ(u_θ，Δt)|，则

θ_k+1＝u_θ

u_θ+1＝θ_k+0.618×(θ_k+1-θ_k)

继续迭代，直到|COV′_θ(l_θ，Δt)|和|COV′_θ(u_θ，Δt)|的差值小于精度ε_θ(一个很小的正数，ε_θ≤10^-6)，即得到θ方向的零点，θ＝(l_θ+u_θ)/2。

(3)计算点(θ，Δt)的COV′，如果COV′的值小于精度ε₃(一个很小的正数，ε₃≤10^-8)，停止循环，将该点作为结果；反之，则继续重复步骤(1)和(2)交替搜索方向，寻找新的Δt和θ方向的零点，直到结果满足条件。

如果一个点在正、负两个方向都可以搜索Δt或θ方向的零点，则需要计算两个方向的零点，并选择较低的COV'值作为更新点。如图4所示，图4a和4b显示了一维进退法求取零点。图4a显示了快速确定零点的范围，图4b显示了通过缩小范围得到零点值。在图4c中，将SWS分析与进退法相结合，寻找真实值。首先求出延时时间方向的零点，然后求出方位角的零点。整个过程搜寻了2次时差方向的零点和2次方位角方向的零点。

本发明将区域划分为四个子区域，并在每个子区域内进行计算，得到四个输出值，比较四个输出值选择最大的COV值对应的点为结果点，以排除局部最优解。根据测试区域为180×60(180为方位角范围，60为时差范围)，所以每个子区域的大小为90×30，共有2821个网格点。图5是进退示意图和区域划分图。取每个子区域的中心点作为初始点，真实值为(30,10)。

根据图6所示，当确定真实值后，可以得到此时的R分量和T分量，然后通过以下公式：

计算得到S1(t)和S2(t)。根据给定的初始点确定区间，通过图3的流程确定该区域的结果值，在其他区域重复此操作，得到4个结果值，比较得到最大的结果值对应的时差和角度作为输出值。

图7是测试范围内每个点依次作为真实值时的合成数据误差图，a为角度误差图，b为时差误差图。

在实际数据测试部分，将B盆地A井电阻率测井得到裂缝孔隙度曲线作为基准线，与在该工区测试的各向异性程度曲线进行对比。A在第50道附近。目标层是S地层(图8，虚线框)。观察到R或T分量中相邻的两个波峰接近32ms。因此，将时间窗口的长度设置为32ms。经过测试的计算统计显示：进退法每次最大迭代次数不超过200次，总共计算次数为25192次，相比旋转相关法所需要计算的3153744次少了2个数量级的计算。本实施例的实际数据计算结果见图9。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims

1.一种基于进退法的横波分裂分析方法，其特征在于，步骤如下：

S1、根据快慢横波相关性构造目标函数如下：

式中，θ和Δt为横波分裂参数，S1(t)和S2(t)分别为一个时窗内的快横波和慢横波的时间序列；和/>分别是S1(t)和S2(t)在一个时窗内的平均值；

S2、通过确定的真实值或实际数据得到R分量数据R(t)和T分量数据T(t)，进一步计算得到S1(t)和S2(t)；

S3、将测试区间划分为4个子区间，取每个子区间的中心点作为初始点；

S4、在每个子区间分别采用二维进退法确定该区间的结果值，共计得到4个结果值；

S5、比较4个结果值，取得到最大COV值对应的时差和角度作为输出值。

2.如权利要求1所述的基于进退法的横波分裂分析方法，其特征在于，所述步骤S3中，总区间为180×60，180为方位角范围，60为时差范围，分为四个子区间后，每个子区域的大小为90×30。

3.如权利要求1所述的基于进退法的横波分裂分析方法，其特征在于，所述步骤S2中，由R分量数据R(t)和T分量数据T(t)，根据以下公式计算得到S1(t)和S2(t)：

4.如权利要求1所述的基于进退法的横波分裂分析方法，其特征在于，所述步骤S4包括以下子步骤：

(1)对Δt方向零点的搜索；

第一步、确定范围，范围的上下界由初始点(θ_k，Δt_k)和比较点(θ_k，Δt_k+1)确定；Δt_k+1等于Δt_k加上步长h；如果COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)小于零，得到范围[Δt_k，Δt_k+1]，进行下一步；如果COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)大于零，得到新的Δt_k+1，新Δt_k+1等于新的Δt_k加上新的h，新的Δt_k等于上一步的Δt_k+1，新的h是上一步h的两倍；然后，再次判断COV′_Δt(θ_k，Δt_k)×COV′_Δt(θ_k，Δt_k+1)的值是大于零还是小于零，确定是执行下一步或者继续更新这两个点；

l_t＝Δt_k+0.382×(Δt_k+1-Δt_k)，

u_t＝Δt_k+0.618×(Δt_k+1-Δt_k).

通过比较|COV′_Δt(θ_k，l_t)|与|COV′_Δt(θ_k，u_t)|的大小，更新l_t或者u_t：

如果|COV′_Δt(θ_k，l_t)|＞|COV′_Δt(θ_k，u_t)|，则

Δt_k＝l_t，

l_t+1＝Δt_k+0.382×(Δt_k+1-Δt_k)

如果|COV′_Δt(θ_k，l_t)|＜|COV′_Δt(θ_k，u_t)|，则

Δt_k+1＝u_t，

u_t+1＝Δt_k+0.618×(Δt_k+1-Δt_k)；

继续迭代，直到|COV′_Δt(θ_k，l_t)|与|COV′_Δt(θ_k，u_t)|之间的差值小于精度ε_Δt，ε_Δt＝10^-6，即得到Δt方向的零点，Δt＝(l_t+u_t)/2；

(2)对θ方向零点的搜索；

第一步、由初始点(θ_k，Δt)得到比较点(θ_k+1，Δt)，θ_k+1等于θ_k加上步长h；然后，计算COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值，如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)＜0，则范围为[θ_k，θ_k+1]，执行下一步；如果COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)＞0，继续更新，得到新的θ_k+1，新的θ_k+1等于新的θ_k加上新的h，新的θ_k，等于上一个θ_k+1，新的h是上一个h的两倍；然后，再次判断COV′_θ(θ_k，Δt)×COV′_θ(θ_k+1，Δt)的值是大于零还是小于零，确定是执行下一步或者继续更新这两个点；

l_θ＝θ_k+0.382×(θ_k+1-θ_k)

u_θ＝θ_k+0.618×(θ_k+1-θ_k)

如果|COV′_θ(l_θ，Δt)|＞|COV′_θ(u_θ，Δt)|，则

θ_k＝l_θ

l_θ+1＝θ_k+0.382×(θ_k+1-θ_k)

如果|COV′_θ(l_θ，Δt)|＜|COV′_θ(u_θ，Δt)|，则

θ_k+1＝u_θ

u_θ+1＝θ_k+0.618×(θ_k+1-θ_k)

5.如权利要求4所述的基于进退法的横波分裂分析方法，其特征在于，所述步骤S4中，如果一个点在正、负两个方向都可以搜索Δt或θ方向的零点，则计算两个方向的零点，并选择较低的COV′值作为更新点。