CN104181590B - 基于小波包分解的叠前道集优化方法 - Google Patents

Abstract

本发明提供一种基于小波包分解的叠前道集优化方法，首先利用小波包将地震道集进行自适应分解，将分解系数划分为信号低频、中频、高频部分，对高频信号进行噪声检测去除，对低频和中频信号分别进行同相轴拉平，最后对拉平后的低频与中频信号进行重构地震道集；本发明的方法将地震信号利用小波包作自适应分解，这种分解无冗余也无疏漏，在各个频段同相轴特征不相同，在每个频段做同相轴相似度计算更为准确；本发明的方法定义吸引度与归属度两类消息进行道集间相似度传播，自动选取参考道和类似道，无需人工干预，可靠性高；本发明的方法以时窗中心点的移动量作为种子点移动量，其他点移动量由种子点移动量插值得到，避免了波形失真。

Description

基于小波包分解的叠前道集优化方法

技术领域

本发明属于地震道集领域，具体涉及一种地震道集拉平方法。

背景技术

页岩气储层裂缝和脆性是影响页岩气产能的最重要的两个因素，储层发育有良好的天然裂缝及层理是实现增产改造的前提条件，脆性指数越高越易形成缝网，有利于增产改造。宽方位角叠前反演可以获得所需的裂缝和脆性属性，但宽方位角叠前反演需要的道集是拉平的，进行道集拉平是反演的前序工作，直接影响反演效果，因此道集拉平具有重要意义。

此外，拉平的道集能够获得好的成像结果，叠前所进行的振幅随偏移距变化(AVO)研究也需要动校正(NMO)之后拉平道集。否则，截距属性和梯度就会受到污染。

目前地震道集的拉平方法主要有：基于速度的拉平方法、基于统计的拉平方法、RSI公司AVATAR软件模块方法等。

基于速度的拉平方法：它假定原始叠前地震道中没有拉平的地震波主要由于RMO(Residual NMO)引起，因此使用二阶或者四阶的RMS速度场的高精度速度估计能够拉平地震波；

基于统计的拉平方法：假定在地震波射线路径上的局部速度场的扰动导致了叠前地震道集的随机起伏，因此不能使用全局的速度场来校正，将它们视作统计误差。

以上两种方法都需要对速度进行准确估计，难度较大，并且上述方法都是基于反射面的测量都是水平的，然而这种假设仅仅适用于合成道集不适用于叠前道集。

RSI公司AVATAR软件模块方法：首先建立L2范数的目标函数，目标函数主要由AVOI或者III类来表示，然后利用时间方向的滑动时窗产生每道的移动解，最后最小化目标函数，对应的解即为最优解；

以上方法仅适用于AVO classI和AVO classIII异常，不适用于AVO classII相位异常，对于AVO classII异常这种方法是不稳定的。基于目前的研究现状本项提出了一种稳定高效的叠前道集优化方法。

Hinkley等人提出了动态道集的拉平算法(DGF)，其思路是通过最小化道集间误差平方和。

该方法是一个最优化处理过程，优点是效率高，适合于处理道集已经基本拉平的数据，不适合运用到同相轴不明显的数据中。

Gulunay等人提出了与动态拉平算法(DGF)类似的道集拉平算法，它是一种简单的道集拉平，这种校正建立在点对点映射的基础上，也就是每个输出数据样点是通过同一道上同样时间上的输入数据样点决定。

该方法实现简单但存在两大缺点：①需要人工选择某道为基准道，然而人工选择基准道是不准确的。②每道只能得到一个整体移动量，整体移动道集后会造成某些同相轴拉平效果较好，某些同相轴反而更差，出现“顾此失彼”现象。

地震道集是通过人工激发地震波，在传播过程中由于波形在地底传播遇到不同介质分界面时发生反射和折射，通过地表的接收机接收返回波形信号形成地震道集，而同相轴是地震记录上各道振动相位相同的极值(俗称波峰或波谷)的连线。拉平的道集能够获得良好的成像效果，是叠前反演工作的前序工作，因此道集分析与优化对于叠前地震数据处理有着非常重要的作用。

发明内容

为解决上述问题，本发明提出了基于小波包分解的叠前道集优化方法，利用小波包对地震道集进行自适应分解，在不同频段进行同相轴拉平后重构地震道集。

本发明的技术方案：基于小波包分解的叠前道集优化方法，具体包括以下步骤：

S1：小波包自适应分解：对地震信号进行小波包分解，去除高频噪声后得到地震低频信号和地震中频信号；

S2：对步骤S1中得到的地震低频信号，进行低频剖面拉平，具体包括以下分步骤：

S21：计算小波系数的欧氏距离：

所述计算小波系数的欧氏距离包括：设叠前道集道数为N，初始化N×N大小的矩阵S，计算同一个时窗内i道与j道低频小波系数的欧氏距离，设Trace1(x₁,x₂,x₃..x_n)为i道数据的低频小波系数，Trace2(y₁,y₂,y₃..y_n)为j道数据的低频小波系数，二者欧式距离计算公式：

$E(i,j)=\sqrt{\underset{k_1=1}{\overset{n}{\Σ}}{(x_{k_1}-y_{k_1})}^2};$

其中，n表示数据点，k₁∈{1,2,3……n}；

S22：记录相似性最大的移动量：

所述记录相似性最大的移动量包括：固定任意一道i的时窗，上下移动道j的时窗，每次移动一个数据点，移动一次记录小波系数的欧式距离，直到找到欧氏距离最小的位置，记录当前移动量，如果是上移，移动量取正，反之，取负；

S23：邻近道传播迭代求取参考道和类似道，具体包括以下分步骤：

S231：定义吸引度矩阵R(i,k)和归属度矩阵A(i,k)，首先初始化为0矩阵，矩阵大小与矩阵S相同，为N×N的矩阵；

S232：更新信息R(i,k)，A(i,k)和A(k,k)；

进一步地，采用更新R(i,k)；

进一步地，采用更新A(i,k)；

进一步地，采用更新A(k,k)；

S233：更新吸引度矩阵R(i,k)和归属度矩阵A(i,k)；

进一步地，采用R(k,k)＝P(k)-max{A(k,i′)+S(k,i′)}更新吸引度矩阵R(i,k)和归属度矩阵A(i,k)；

S234：确定道集i₁的类似道k；

S235：确定参考道i₂与类似道k后，选择C矩阵的第i行，C(i,i)作为当前时窗参考道i₂的移动量，C(i,j)作为当前时窗类似道的移动量，其余地震道根据参考道与类似道的移动量插值得到；

S24：求每个点移动量：根据每个时窗中心点的移动量线性插值得到每个点的移动量，实施道集优化；

S3：对步骤S1中得到的中频信号，进行中频剖面拉平；

S4：重构地震信号，根据拉平后的低频和中频小波包系数重构地震信号；

更进一步地，更新后的吸引度矩阵R(i,k)和归属度矩阵A(i,k)为：

$\left\{\begin{array}{c}{R}_{new}(i,k)=\mathrm{\λ}\×R_{old}(i,k)+(1-\mathrm{\λ})\×R(i,k)\\ A_{new}(i,k)=\mathrm{\λ}\×A_{old}(i,k)+(1-\mathrm{\λ})\×A(i,k)\end{array}\right..$

本发明的基于小波包分解的叠前道集优化方法，利用小波包分解地震道集，在各个频段上，最大化道集间的相似性确定时窗中心点移动量，在道集间传递吸引度与归属度两类消息实现参考道和类似道的选取，具有如下优点：

1)将地震信号利用小波包作自适应分解，这种分解无冗余也无疏漏，在各个频段同相轴特征不相同，在每个频段做同相轴相似度计算更为准确；

2)定义吸引度与归属度两类消息进行道集间相似度传播，自动选取参考道和类似道，无需人工干预，可靠性高；

3)以时窗中心点的移动量作为种子点移动量，其他点移动量由种子点移动量插值得到，避免了波形失真。

附图说明

图1叠前道集分析与优化流程图。

图2一维三阶小波包分解流程图。

图3一维三阶小波包分解频率特征。

图4原始地震信号。

图5小波包分解后的高频、中频、低频分解系数。

图6原始地震信号、去噪后地震信号与噪声示意图。

图7地震剖面低频信号与中频信号示意图。

图8低频地震剖面拉平算法流程。

图9根据欧式距离寻找最优移动量。

图10吸引度与归属度传递过程。

图11低频剖面拉平前与拉平后示意图。

图12中频地震剖面拉平算法流程。

图13中频地震剖面拉平前与拉平后示意图。

图14原始数据与道集优化结果示意图。

具体实施方式

下面结合附图对本发明的基于小波包分解的叠前道集优化方法，做进一步的说明。

本发明首先利用小波包将地震道集进行自适应分解，将分解系数划分为信号低频、中频、高频部分，对高频信号进行噪声检测去除，对低频和中频信号分别进行同相轴拉平，最后对拉平后的低频与中频信号进行重构地震道集。总流程图如图1所示，具体方案如下

S1：小波包自适应分解：对地震信号进行小波包分解得的地震低频信号和地震中频信号；

小波包变换的意义就在于能够在不同尺度上对信号进行分解，图2与图3分别是一维三阶小波包分解流程图及频率特征图，相对小波分解，小波包分解能进一步细分高频部分，达到更优的时频局部化效果。提取一道地震数据进行分解，其分解示意图如图4。本发明对道集做两层分解得到低频、中频、高频小波分解系数，如图5。

去除高频噪声后得到去噪前后的剖面，如图6，左图为原始地震信号剖面，中图为去噪后地震信号剖面，右图为噪声，由噪声剖面可以看出，去噪后的剖面同相轴更为清晰，同相轴结构特征保持较好。

地震信号的低频和中频成分，如图7，左图为地震低频信号，右图为地震中频信号。

S2：低频剖面拉平：地震信号低频部分包含了地震剖面主要的同相轴，设置参数窗口大小N_w，窗口移动量N_s，通过计算地震剖面的主频，定义窗口大小为2个信号周期，窗口移动量为一个信号周期。固定任意一道时窗，上下移动其他道窗口，求取相似度最大时的移动量，通过近邻相似传播迭代求取参考道与相似道，最后根据参考道和相似道移动量求取剩余道的移动量。其算法流程如图8所示：

S21：计算小波系数的欧氏距离：设叠前道集道数为N，初始化N×N大小的矩阵S，计算同一个时窗内i道与j道低频小波系数的欧氏距离，设Trace1(x₁,x₂,x₃..x_n)为i道的低频小波系数，Trace2(y₁,y₂,y₃..y_n)为j道数据的低频小波系数，二者欧式距离计算公式：

$E(i,j)=\sqrt{\underset{k_1=1}{\overset{n}{\Σ}}{(x_{k_1}-y_{k_1})}^2};$

其中，n表示数据点，k₁∈{1,2,3……n}；

S22：记录相似性最大的移动量：固定任意一道i的时窗，上下移动道j的时窗，每次移动一个数据点，移动一次记录小波系数的欧式距离，直到找到欧氏距离最小的位置，记录当前移动量，如果是上移，移动量取正，反之，取负。图9显示了两道数据同一时窗内移动量与欧氏距离的关系，相似性最大点移动量为9(即向下移动9个数据点)。初始化N×N大小的矩阵C，矩阵C(i,j)每一个行向量表示固定i道时窗，移动j道数据时窗对应的最优移动量。矩阵S记录相应的最小欧式距离。

S23：邻近道传播迭代求取参考道和类似道如图10所示：叠前道集优化需要确定参考道和类似道，类似道是是指同相轴与参考道非常近似的地震道，在确定了参考道和类似道后，剩余的地震道移动量可以根据参考道和类似道求得，因此参考道与类似道的求取变得尤为重要。本发明通过欧氏距离矩阵S作为原始数据，定义吸引度矩阵R(i,k)和归属度矩阵A(i,k)表示道集之间的两类信息，其中R(i,k)是从道集i₁指向候选类似道k，它反映了道集k适合作为道集i₁的参考道集的合适程度；A(i,k)从候选类似道k指向道集i₁，它反映了道集i₁选择k作为其参考道的合适程度，在两类消息传播结束时，道集i₁的参考道集确定为k，k满足：

$\underset{k}{\mathrm{argmax}}\left(A\right(i,k)+R(i,k\left)\right)$

为了便于计算将S矩阵中的每个最小欧氏距离取负号，定义吸引度矩阵R(i,k)和归属度矩阵A(i,k)两个目标函数，首先初始化为0矩阵，矩阵大小与矩阵S相同，为N×N的矩阵，迭代求解R(i,k)和A(i,k)：

根据更新信息R(i,k)。

根据更新信息A(i,k)。

根据更新信息A(k,k)。

根据式R(k,k)＝P(k)-max{A(k,i′)+S(k,i′)}更新吸引度与归属度：

$\left\{\begin{array}{c}{R}_{new}(i,k)=\mathrm{\λ}\×R_{old}(i,k)+(1-\mathrm{\λ})\×R(i,k)\\ A_{new}(i,k)=\mathrm{\λ}\×A_{old}(i,k)+(1-\mathrm{\λ})\×A(i,k)\end{array}\right.$

其中，P(k)取S矩阵的均值，λ是收敛系数，主要用于调节算法的收敛速度与迭代过程的稳定性，λ(0<λ<1)为收敛系数，λ越大消除振荡的效果越好，但收敛速度也越慢，反之亦然。下标old和new分别代表上一次和本次更新消息的最终结果。迭代终止条件为迭代次数大于预先设定次数或者迭代数次R(i,k)与A(i,k)保持不变。

确定道集i的类似道k：

$k=\arg \underset{k}{max}\{A(i,k)+R(i,k)\}$

运用上式，若i＝k,则道集i是参考道集；若i≠k,则道集k是道集i的参考道集，以此将道集分成几组，将每组中地震道对应的相似度矩阵C所在行求和，取和值最小的一组即相似度最大的一组，没有选中的道集作为剩余道。剩余道的移动量根据参考道和类似道的移动量插值得到。从而得到每个时窗中心点的移动量。滑动时窗，重复上面的步骤S21、步骤S22、步骤S23，得到每个时窗中心点的移动量。图8是参考道与类似道选择流程。

确定参考道i与类似道k后，选择C矩阵的第i行，C(i,i)作为当前时窗参考道i的移动量，C(i,j)作为当前时窗类似道的移动量，其余地震道根据参考道与类似道的移动量插值得到。例如：第5道为参考道，其移动量为-4，第8道为类似道，其移动量为+3，中间的6、7道为剩余道，则6、7道的移动量为-5/3，7/3。

S24：求每个点移动量：根据每个时窗中心点的移动量线性插值得到每个点的移动量，实施道集优化。图11为截取低频剖面一部分显示拉平前后的效果，从效果上看低频剖面同相轴拉平效果很好，同相轴得到很好拉平。

S3:中频剖面拉平：地震道集中频剖面拉平算法与低频拉平算法思想类似，相对于低频同相轴清晰且“粗壮”，中频信号同相轴比较密集，分辨率较低。中频剖面拉平算法需要修改两点：(一)窗口大小与窗口移动量需要重新设置，具体设置方式根据中频剖面的主频确定，中频剖面主频比低频剖面主频高，因此一个信号周期中点数相较于低频剖面更少，因此窗口大小及窗口移动量适当调小(二)因中频同相轴比较密集，在求取各道时窗中心点的移动量时尽可能多的选择参考道，类似道，减少插值道数，防止信号失真。因此修改目标函数R(k,k)中的P(k)，取S矩阵最小值，迭代得到组数最少，选择组内道数最多的一组。中频剖面拉平算法如图12，图13为截取中频剖面一部分显示拉平前后的效果，从效果上看拉平前中频剖面同相轴比较乱，通过自动拉平，剖面主要同相轴比较平整。

S4：重构地震信号：根据拉平后的低频和中频小波包系数重构地震信号。图14是道集优化的结果。从结果上看通过小波包对地震信号进行自适应分解，针对不同频段同相轴的特点对各频段小波系数进行不同处理后重构地震信号，这套方案能有效地进行道集优化，同时完成了地震数据的去噪及同相轴的拉平，是一种稳定高效的道集优化方法。

本发明利用小波包分解地震道集，在各个频段上，最大化道集间的相似性确定时窗中心点移动量，在道集间传递吸引度与归属度两类消息实现参考道和类似道的选取，提出了一种稳定高效的叠前道集优化方法，具有如下优点：

(1)将地震信号利用小波包作自适应分解，这种分解无冗余也无疏漏，在各个频段同相轴特征不相同，在每个频段做同相轴相似度计算更为准确。

(2)定义吸引度与归属度两类消息进行道集间相似度传播，自动选取参考道和类似道，无需人工干预，可靠性高。

(3)以时窗中心点的移动量作为种子点移动量，其他点移动量由种子点移动量插值得到，避免了波形失真。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (3)

1.基于小波包分解的叠前道集优化方法，其特征在于，包括：

S1：小波包自适应分解：对地震信号进行小波包分解，去除高频噪声后得到地震低频信号和地震中频信号；

S2：对步骤S1中得到的地震低频信号，进行低频剖面拉平，具体包括以下分步骤：

S21：计算小波系数的欧氏距离：

所述计算小波系数的欧氏距离包括：设叠前道集道数为N，初始化N×N大小的矩阵S，计算同一个时窗内i道与j道低频小波系数的欧氏距离，设Trace₁(x₁,x₂,x₃..x_n)为i道数据的低频小波系数，Trace₂(y₁,y₂,y₃..y_n)为j道数据的低频小波系数，二者欧式距离计算公式：

$E(i,j)=\sqrt{\underset{k_1=1}{\overset{n}{\&Sigma;}}{(x_{k_1}-y_{k_1})}^2};$

其中，n表示数据点，k₁∈{1,2,3……n}；

S22：记录相似性最大的移动量：

所述记录相似性最大的移动量包括：固定任意一道i的时窗，上下移动道j的时窗，每次移动一个数据点，移动一次记录小波系数的欧式距离，直到找到欧氏距离最小的位置，记录当前移动量，如果是上移，移动量取正，反之，取负；

S23：邻近道传播迭代求取参考道和类似道，具体包括以下分步骤：

S231：定义吸引度矩阵R(i,k)和归属度矩阵A(i,k)，首先初始化为0矩阵，矩阵大小与矩阵S相同，为N×N的矩阵；

S232：更新信息R(i,k)，A(i,k)和A(k,k)；

采用更新R(i,k)；

采用更新A(i,k)；

采用更新A(k,k)；

S233：更新吸引度矩阵R(i,k)和归属度矩阵A(i,k)；

采用R(k,k)＝P(k)-max{A(k,i′)+S(k,i′)}更新吸引度矩阵R(i,k)和归属度矩阵A(i,k)；

其中，P(k)取S矩阵的均值；

S234：确定道集i₁的类似道k；

S235：初始化N×N大小的矩阵C，矩阵C(i,j)每一个行向量表示固定i道时窗，移动j道数据时窗对应的最优移动量；

确定参考道i₂与类似道k后，选择C矩阵的第i行，C(i,i)作为当前时窗参考道i₂的移动量，C(i,j)作为当前时窗类似道的移动量，其余地震道根据参考道与类似道的移动量插值得到；

S24：求每个点移动量：根据每个时窗中心点的移动量线性插值得到每个点的移动量，实施道集优化；

S3：对步骤S1中得到的中频信号，进行中频剖面拉平；

S4：重构地震信号，根据拉平后的低频和中频小波包系数重构地震信号。

2.根据权利要求1基于小波包分解的叠前道集优化方法，其特征在于，更新后的吸引度矩阵R(i,k)和归属度矩阵A(i,k)为：

$\left\{\begin{array}{c}{R}_{new}(i,k)=\mathrm{\&lambda;}\&times;R_{old}(i,k)+(1-\mathrm{\&lambda;})\&times;R(i,k)\\ A_{new}(i,k)=\mathrm{\&lambda;}\&times;A_{old}(i,k)+(1-\mathrm{\&lambda;})\&times;A(i,k)\end{array}\right.;$

其中，λ是收敛系数，R_new(i,k)表示本次迭代的吸引度矩阵、A_new(i,k)表示本次迭代的归属度矩阵，R_old(i,k)表示上一次迭代的吸引度矩阵，A_old(i,k)表示上一次迭代的归属度矩阵；R(i,k)表示更新后的吸引度矩阵，A(i,k)表示更新后的归属度矩阵。

3.根据权利要求1基于小波包分解的叠前道集优化方法，其特征在于，所述S234确定道集i₁的类似道k的公式为：

$k=\arg \underset{k}{max}\{A(i,k)+R(i,k)\};$

若i＝k,则道集i₁是参考道集；若i≠k,则k是道集i₁的类似道。