CN104570074A - 一种基于奇异值分解技术的废道识别方法 - Google Patents

Abstract

本发明提供了一种基于奇异值分解技术的废道识别方法，属于石油地球物理勘探领域。本方法包括：(1)将地震数据剖面划分为三个范围，分别为浅范围、中范围和深范围，其中浅范围包括初至前的数据，中范围包括初至所在区域内的数据，深范围包括初至后的深层数据；(2)对于每个地震道的数据，分别获取其三个范围中的属性特征，所述属性特征包括浅范围的平均振幅、主频、过零点个数，中范围的平均振幅、主频、过零点个数、相关系数，深范围的平均振幅、主频、过零点个数、视衰减因子；(3)对每个属性特征分别作归一化处理，得到属性矩阵；(4)对所述属性矩阵进行奇异值分解，并计算得到欧氏距离；(5)输入阈值。

Description

一种基于奇异值分解技术的废道识别方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种基于奇异值分解技术的废道识别方法。

背景技术

在地震资料处理过程中，废道检测与剔除是数据处理中非常重要且繁琐的工作。通常处理人员只能通过振幅值的大小、视频率的高低等特征识别废道，人工判断和手工剔除存在着效率低、不够准确的问题，难以适应高密度勘探中海量地震数据处理进度的要求。

目前一些自动或半自动地震道编辑算法和软件，例如软件Promax，它的处理过程主要是机械地使用该软件中的各单独处理模块，通过主频来判别异常道，因此处理结果通常不尽如人意，最终仍需与人工编辑相结合。有的采用分时间段数据分析，提取不同时间段内地震数据的特征参数，综合互相关参数、分频扫描、视衰减因素等检测不正常道。然而分析过程中根据各属性特点的不同设定不同的判定条件，算法比较复杂，不容易实现。而一些基于地震道特征评价的识别方法，由于其参数众多或准则复杂，不利于大规模推广应用。

考虑到提取的地震道特征属性是用矩阵来描述的，本发明将奇异值分解算法引入废道识别算法中，由于矩阵的奇异值是矩阵特征值的改进，用奇异值代表矩阵的特征比用特征值更有优点，可以应用于一般矩阵的求解，而不仅仅是方阵。同时降低存储矩阵的大小，有利于减少内存空间的占用。

在现今的海量地震数据处理中，废道识别依然成为地震资料处理中的一个瓶颈。由于废道对后续处理效果的影响不容忽视，特别是能量较强的废道，对地表一致性处理振幅补偿、叠前随机噪声衰减等多道处理有不良的影响：影响振幅的衰减规律，导致补偿系数难以准确求取；使得地震反射信号在频率域的可预测性减弱，限制了叠前随机噪声衰减作用的有效发挥。

在现有的废道识别技术中，0mega作为目前主流的处理软件，在滤波过程中可以识别出工业干扰、空炮等废道，但仍旧无法识别出夹杂在正常数据中的空道；在Promax中，只能采取人工操作的方式来识别废道，费时费力。然而由于废道对地震资料处理的结果有着很大的影响，对这些废道的处理也势在必行，而仅靠人工拾取废道远不能满足现有的项目进度要求。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于奇异值分解技术的废道识别方法，高效使用的地震道自动编辑算法可以进一步减轻处理人员的工作量，并且可以更准确的进行废道的剔除，便于后续处理工作的进行。

本发明是通过以下技术方案实现的：

一种基于奇异值分解技术的废道识别方法，包括：

(1)将地震数据剖面划分为三个范围，分别为浅范围、中范围和深范围，其中浅范围包括初至前的数据，中范围包括初至所在区域内的数据，深范围包括初至后的深层数据；

(2)对于每个地震道的数据，分别获取其三个范围中的属性特征，所述属性特征包括浅范围的平均振幅、主频、过零点个数，中范围的平均振幅、主频、过零点个数、相关系数，深范围的平均振幅、主频、过零点个数、视衰减因子；

(3)对每个属性特征分别作归一化处理，得到属性矩阵；

(4)对所述属性矩阵进行奇异值分解，并计算得到欧氏距离；

(5)输入阈值；

(6)判断所述欧氏距离是否大于阈值，如果是，则此地震道的识别结果为废道，如果否，则此地震道的识别结果为正常道。

所述方法进一步包括：

(7)判断所述识别结果是否准确，如果否，则调整阈值，然后返回步骤(6)；如果是，则转入步骤(8)；

(8)结束。

所述步骤(1)中的所述浅范围和深范围均为矩形区域，所述中范围为倒V型区域。

所述步骤(2)中获取其三个范围中的地震属性是这样实现的：

提取浅范围内的地震道的平均振幅、主频和过零点个数；

提取中范围内的地震道的平均振幅、主频和过零点个数，并计算得到中范围内相邻道的相关系数；

提取深范围内的地震道的平均振幅、主频、过零点个数，并计算得到该范围内的视衰减因子；

这样得到11个地震属性。

所述步骤(3)是这样实现的：

对于含有n道数据的炮集，对于第j道数据，将所述11个属性特征，记为a_i＝[a_1j，a_2j，a_3j，a_4j，a_5j，a_6j，a_7j，a_8j，a_9j，a_10j，a_11j]^T，其中a_1j，a_2j，a_3j表示浅范围内的平均振幅、主频和过零点个数；a_4j，a_5j，a_6j，a_7j表示中范围内的平均振幅、主频、过零点个数和相关系数；a_8j，a_9j，a_10j，a_11j表示深范围内的平均振幅、主频、过零点个数和视衰减因子；

对每个属性特征分别作归一化处理，得到属性矩阵

所述步骤(4)是这样实现的：

对所述属性矩阵A采用奇异值分解，得到右特征矩阵V，取右特征矩阵V的第一列v₁，即v₁＝[v₁₁，…，v_n1]^T，得到列向量s＝Av₁，即

利用下式求得属性矩阵A中的每一列向量与列向量s的欧氏距离d_j：

$d_j=\left|a_j\right|\×\sqrt{1-{\left(\frac{a_j\·s}{\left|a_j\right|\×\left|s\right|}\right)}^2}.$

所述阈值的取值范围为[0，1]之间。

所述步骤(7)中的调整阈值是这样实现的，

如果地震道中的废道没有被识别出，则将阈值减小，如果将正常道判定为了废道，则将阈值增大。

与现有技术相比，本发明的有益效果是：本发明方法实现了对相关废道的快速识别，可以自动的拾取废道(包含空道，噪声道，野值道等)，拾取结果准确，速度较快；此外本发明可以作为一个单独的自动化模块进行废道剔除，从而推进整个数据处理的进度和精度。

附图说明

图1是根据不同时间段确定不同识别废道区域的示意图；

图2是本发明实施例中的炮号1的识别废道效果示意图；

图3是图2中部分剖面的放大图；

图4是本发明实施例中的炮号195的识别废道效果示意图；

图5是图4中部分剖面的放大图；

图6是本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

本发明从废道与正常道的相异处入手，提取出包括平均振幅、主频、过零点个数、相关系数、视衰减因子等11项特征，借助奇异值分解算法识别异常特征，进而识别出废道，本发明简单易懂，内存使用量较少，可以实现海量地震数据的废道识别。

(1)原始地震道属性特征提取

初至前数据的分析与检测。地震波初至到达前，检波器接收到的是环境噪声的信号，因此初至前的地震数据可用于环境噪声的分析和不正常道的识别。通常情况下，假设环境噪声为平稳随机过程，可通过对初至前的数据进行统计来研究环境噪声的规律。对于单炮数据，此处选择初至前的一个矩形区域，选择平均振幅、主频和过零点个数三个特征作为识别废道的属性特征。图1中的数据采样间隔为2ms，取前200个采样点进行初至前数据的分析。

在初至所在区域内的废道特征提取。对于原始单炮数据，如果划定初至所在的一个较大范围，可以得到一个“倒V”型的区域，相对于原始单炮数据上的“矩形”区域，它可以尽量避免由于选择区间不当造成的属性特征的不同。在此区域中，提取地震道的平均振幅、主频、过零点个数三个特征；此外，在此区域中含有品质较好的反射层，可以计算此区域中相邻道的相关系数(依据地震道的相关系数(或者称为互相关函数)公式得到，具体步骤为：1)找到零偏移距对应的地震信号x_k(n)，0≤n≤M；2)对于地震信号x_i(n)，0≤n≤L，查看其偏移距offset_i，若offset_i≥0，则第x_i(n)道与x_k(n)的互相关函数为r_ik(n)，否则第x_i(n)道与x_k(n)的互相关函数为r_ki(n)；3)无论r_ik(n)还是r_ki(n)都是利用FFT(快速傅里叶变换)计算得到的；具体为选择N＝2^k，N≥L+M，构造x_k(n)，x_i(n)的以N为周期的周期信号利用FFT分别计算对应的离散频谱计算或者在利用FFT反变换得到r_ik(n)或r_ki(n)；4)在0≤n≤N-L寻找r_ik(n)或r_ki(n)的最大值，即为第i道的相关系数。)，如果是正常道，数据有良好的相关性，反之，数据中有异常道，其数据的相关性会发生变化，其相关最大值及其对应的时移量都会出现异常。因此，在此区间选择上述四项属性特征识别废道。图1中的数据取600个采样点进行初至区域的数据分析。

深层区域内的废道特征提取。当地震波信号传播到深层后，由于大地滤波的作用，地震波的能量、振幅等都会呈现指数衰减，然而，大部分的废道却不会出现能量的衰减，所以，在深层的矩形区域内，选择平均振幅、主频、过零点个数以及视衰减因子(此处计算的是视衰减因子，具体为：对于第i道地震数据，其视衰减因子为其中A(t₁)为t₁时刻的振幅，在本发明中，A(t₁)为深层起始位置的振幅的平方，根据用户的交互操作得到，A(t₂)为深层结束位置的振幅的平方)作为废道识别的属性特征。图1中的数据取最后400个采样点进行深层区域的分析。

具体实施时，这三个范围是用户通过交互的方式在界面上获取的。

(2)采用奇异值分解算法识别废道

奇异值分解是线性代数中一种重要的矩阵分解，是矩阵分析中正规矩阵酉对角化的推广，在信号处理、统计学等领域有重要作用。它在统计中的主要应用为主成分分析，它是一种数据分析方法，用来找到大量数据中所隐含的“模式”(即规律特点)，起到数据降维的作用。在废道识别中，可以发现正常道具有相似的“模式”，而废道具有相异的“模式”，因此可以利用奇异值分解技术对地震道的不同“模式”进行识别。

引理(奇异值分解)

若矩阵A∈R^m×n，则存在正交矩阵U＝[u₁,u₂，…，u_m]∈R^m×m，V＝[v₁，v₂，…，v_m]∈R^v×v使得U^TAV＝diag[σ₁，σ₁，…，σ_p]＝∑，p＝min(m,n)，则称A＝U∑V^T为矩阵A的奇异值分解。其中σ₁≥σ₂≥…≥σ_p≥0，σ_i(i＝1，2，…，p)为A的奇异值。

通过引理可以发现，奇异值分解是一个适用于任意矩阵的分解方法，矩阵U是矩阵AA^T的特征向量，矩阵V是矩阵A^TA的特征向量，矩阵AA^T和A^TA都有非常清楚的物理意义。由于矩阵AA^T中的每一个元素表示对应地震属性特征的相关性，数值越大相关性越大，矩阵A^TA中的每一个元素表示对应地震道之间的相关性，数值越大相关性越大。在本发明中，根据矩阵与地震道属性的关系，借助奇异值分解技术识别废道。

(3)废道识别算法

如图6所示，具体包括：

第一步：根据地震道属性特征简介，对于含有n道数据的炮集，对于第j道数据，分别提取不同时间段的共11个属性特征，记为a_i＝[a_1j，a_2j，a_3j，a_4j，a_5j，a_6j，a_7j，a_8j，a_9j，a_10j，a_11j]^T，其中a_1j，a_2j，a_3j表示初至前的平均振幅、主频、过零点个数；a_4j，a_5j，a_6j，a_7j表示中间“倒V”型区域中的平均振幅、主频、过零点个数，相关系数；a_8j，a_9j，a_10j，a_11j表示深层处的平均振幅、主频、过零点个数，视衰减因子；

图1中需要的只是三个范围，起始位置到第一条浅灰色线条为第一个范围即浅层范围，第二个范围为初至附近的“倒V”型区域，第三个范围为最后一条浅灰色线条到最后即深层范围。

具体实施时，用户是这样得到这三个范围的：

划分这三个范围时需要用户交互得到，具体步骤：为得到第一个范围，用户需在地震数据剖面上初至上方的某一个位置点击一下，界面上出现一条横线以及一个对话框，“是否确定第一个范围”，点击“确定”，划分出第一个范围；为得到第二个范围，用户在初至附近(初至的上侧或下恻)点击得到两个点，界面上会出现“倒V”型区域，以及一个对话框，“是否确定第二个范围”，点击确定，得到第二个范围；为得到第三个范围，用户在深层区域某一位置点击一下，界面上出现一条横线以及一个对话框，“是否确定第三个范围”，点击确定得到第三个范围。

第二步：将n道数据的属性特征提出后，对每个属性特征分别作归一化处理(即将最小的看做为0，最大的看做为1，其余的按比例得到位于[0，1]之间的数值)，得到最后的矩阵

第三步：对矩阵A采用奇异值分解，得到右特征矩阵V(V是一个n×n的方阵，这是奇异值分解过程中用到的矩阵，具体表达式为取右特征矩阵V的第一列v₁，即v₁＝[v₁₁，…，v_n1]^T，得到列向量s＝Av₁，即

第四步：求矩阵A中的每一列向量与列向量s的欧氏距离d_j，设定阈值d₀，若d_i＜d₀，则j为正常道，否则，为废道。

其中 $d_j=\left|a_j\right|\×\sqrt{1-{\left(\frac{a_j\·s}{\left|a_j\right|\×\left|s\right|}\right)}^2}.$

式中的点表示两个向量之间的乘积，叉表示两个实数之间的乘积；aj指的是第j个列向量，ai表示第i个列向量，aij表示的是矩阵第i行，第j列位置的元素，是一个实数。

考虑到不同信噪比的资料对应不同的阈值，为了使本方法的结果更完美，交互操作过程还包括以下几步：

第一步，查看需要处理的资料，根据信噪比的高低给出一个阈值，它在[0，1]之间，对于信噪比高的资料，对应的阈值较高，信噪比较低的资料，对应的阈值较低。

第二步，根据识别的结果(用户从界面上得到的显示效果，其中废道会用红色标识出来。)调整阈值的大小(如果地震道中的废道没有被识别出，将阈值调小，如果地震道中将正常道视为废道，则将阈值增大)，最后得到合适的识别效果。

本发明中将奇异值分解技术引入到废道识别中，下面以某些地区的实例来说明发明效果。首先给出单炮数据中含异常道较少的地震道，根据本发明设置阈值为0.76，其中图2中含有强噪声道第198道，空道第209道和214道，可以看出本方法可以检测出这三道，图3为图2部分剖面的放大图。然后给出单炮数据中含异常道较多的地震炮集数据，设定阈值为0.7，图5中第238-275，470-500道均为空道，第25道和196道为强噪声道，从图中显示可以看出，这些道都可以识别出来，图5为图4部分剖面的放大图。

本发明针对废道的属性特点，分不同时间段提取包括平均振幅、主频、过零点个数、相关系数、视衰减因子等共11个属性特征，得到地震道属性特征矩阵，然后利用奇异值分解技术，通过设定阈值识别废道。紧接着给出实际资料进行检验，可以看出本算法可以有效的识别出废道，由于奇异值技术的引用，本发明时间复杂度也相应降低。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (7)

1.一种基于奇异值分解技术的废道识别方法，其特征在于：所述方法包括：

(1)将地震数据剖面划分为三个范围，分别为浅范围、中范围和深范围，其中浅范围包括初至前的数据，中范围包括初至所在区域内的数据，深范围包括初至后的深层数据；

(2)对于每个地震道的数据，分别获取其三个范围中的属性特征，所述属性特征包括浅范围的平均振幅、主频、过零点个数，中范围的平均振幅、主频、过零点个数、相关系数，深范围的平均振幅、主频、过零点个数、视衰减因子；

(3)对每个属性特征分别作归一化处理，得到属性矩阵；

(4)对所述属性矩阵进行奇异值分解，并计算得到欧氏距离；

(5)输入阈值；

(6)判断所述欧氏距离是否大于阈值，如果是，则此地震道的识别结果为废道，如果否，则此地震道的识别结果为正常道。

2.根据权利要求1所述的基于奇异值分解技术的废道识别方法，所述方法进一步包括：

(7)判断所述识别结果是否准确，如果否，则调整阈值，然后返回步骤(6)；如果是，则转入步骤(8)；

(8)结束。

3.根据权利要求2所述的基于奇异值分解技术的废道识别方法，其特征在于：所述步骤(1)中的所述浅范围和深范围均为矩形区域，所述中范围为倒V型区域。

4.根据权利要求3所述的基于奇异值分解技术的废道识别方法，其特征在于：所述步骤(2)中获取其三个范围中的地震属性是这样实现的：

提取浅范围内的地震道的平均振幅、主频和过零点个数；

提取中范围内的地震道的平均振幅、主频和过零点个数，并计算得到中范围内相邻道的相关系数；

提取深范围内的地震道的平均振幅、主频、过零点个数，并计算得到该范围内的视衰减因子；

这样得到11个地震属性。

5.根据权利要求4所述的基于奇异值分解技术的废道识别方法，其特征在于：

所述步骤(3)是这样实现的：

对于含有n道数据的炮集，对于第j道数据，将所述11个属性特征，记为a_i＝[a_1j，a_2j，a_3j，a_4j，a_5j，a_6j，a_7j，a_8j，a_9j，a_10j，a_11j]^T，其中a_1j，a_2j，a_3j表示浅范围内的平均振幅、主频和过零点个数；a_4j，a_5j，a_6j，a_7j表示中范围内的平均振幅、主频、过零点个数和相关系数；a_8j，a_9j，a_10j，a_11j表示深范围内的平均振幅、主频、过零点个数和视衰减因子；

对每个属性特征分别作归一化处理，得到属性矩阵

6.根据权利要求5所述的基于奇异值分解技术的废道识别方法，其特征在于：

所述步骤(4)是这样实现的：

对所述属性矩阵A采用奇异值分解，得到右特征矩阵V，取右特征矩阵V的第一列v₁，即v₁＝[v₁₁，…，v_n1]^T，得到列向量s＝Av₁，即

利用下式求得属性矩阵A中的每一列向量与列向量s的欧氏距离d_j：

$d_j=\left|a_j\right|\×\sqrt{1-{\left(\frac{a_j\·s}{\left|a_j\right|\×\left|s\right|}\right)}^2}.$

所述阈值的取值范围为[0，1]之间。

7.根据权利要求6所述的基于奇异值分解技术的废道识别方法，其特征在于：所述步骤(7)中的调整阈值是这样实现的，

如果地震道中的废道没有被识别出，则将阈值减小，如果将正常道判定为了废道，则将阈值增大。