基于波原子分解自适应滤波的地震勘探信号增强方法与流程

[0001]
本发明属于地震勘探技术领域，具体涉及一种基于波原子分解自适应滤波的地震勘探信号增强方法。


背景技术：

[0002]
地震波在地下介质的传播过程中会产生衰减，实际采集到的地震数据中，有效信号非常微弱，加之受到各种噪声的影响使其更加微弱，甚至被打断、畸变或者湮没。随机噪声作为背景噪声，对有效信号的识别与提取影响较大，因此，相应的滤波方法应运而生并不断得到改进，性能优越的方法将会使随机噪声大大消减而有效信号显著增强，从而更有利于后续的资料解释与分析工作。
[0003]
当前地震勘探资料处理中有众多的滤波技术，主要包括：
[0004]
1、通过估计信号特征设计的信号增强方法，如f-x预测滤波；
[0005]
2、利用相关性进行有效信号增强的滤波方法，如k-l变换法；
[0006]
3、基于二项式时频变换的滤波方法，如基于wigner-ville分布(wvd)和pseudo wigner-ville分布(pwvd)的时频峰值滤波；
[0007]
4、基于多尺度分解的滤波方法，如小波滤波、曲波滤波等。
[0008]
虽然这些方法都在随机噪声消减以及有效信号增强方面取得了显著的效果，但是也存在固有的缺点使得滤波效果的提升比较有限。
[0009]
那么如何更好地从随机噪声环境中提取地震勘探有效信号，是目前亟待解决的问题。


技术实现要素：

[0010]
本发明的核心思路为：由于波原子变换具有很好的方向性，且对纹理图像具有最优的稀疏表示，因此，利用波原子变换对含噪的纹理图像进行分解，能够很好地得到纹理方向上的特征信息，再通过自适应维纳滤波(wiener)对分解所得系数进行处理，即可滤除噪声系数而保留有效信号系数，最后通过波原子重构得到增强后的纹理图像。
[0011]
本发明旨在提供一种基于波原子分解自适应滤波的地震勘探信号增强方法，其对经过循环平移的含噪地震数据进行波原子分解后，对所得分解系数进行自适应wiener滤波，然后对滤波后的系数进行波原子重构，最后经过反循环平移得到有效信号增强的地震数据。
[0012]
为实现本发明目的，本发明的技术解决方案为：
[0013]
基于波原子分解自适应滤波的地震勘探信号增强方法，其特征在于，包括以下步骤：
[0014]
步骤1：选取一幅二维含噪地震记录，对其数据进行二维循环平移；
[0015]
步骤2：对循环平移后的含噪记录进行二维波原子分解，得到波原子分解系数；
[0016]
步骤3：对分解所得系数进行自适应wiener滤波，得到滤波后的波原子系数；
[0017]
步骤4：对滤波后的波原子系数进行重构，再进行反循环平移，最终得到有效同相轴增强的地震记录。
[0018]
进一步地，步骤1的具体操作步骤包括：
[0019]
s11：定义选取的含噪地震记录为u(x1,x2)；
[0020]
s12：对含噪地震记录数据进行二维循环平移，且平移依据的公式为：
[0021][0022]
其中，为平移算子，h1，h2为网格上的平移量；
[0023]
对数据进行二维循环平移能够改善在图像边缘或纹理等不连续点邻域内产生伪吉布斯现象而导致有效信号失真的情况；
[0024]
进一步地，步骤2的具体操作步骤包括：
[0025]
s21：定义波原子为其中μ＝(j,m,n)，m＝(m1,m2)和n＝(n1,n2)，并且参量j,m1,m2,n1,n2为整数；
[0026]
s22：定义相空间中一点(x
μ
,ω
μ
)，其中：位置向量x
μ
＝2-j
n，表征的空域中心，波向量ω
μ
＝π2
j
m，表征的频域中心，且c1、c2为两个正常数，则基于波包框架的波原子为：
[0027][0028][0029]
s23：定义实值波包将波原子系数c
j,m,n
作为尺度2-j
上的抽样卷积，可得：
[0030][0031]
记hilbert变换为h，定义正交基及对偶正交基分别为：
[0032][0033][0034]
记及则形成了波原子紧支框架，且有：
[0035][0036]
根据上式可得，二维波原子分解系数为：
[0037][0038]
进一步地，步骤3的具体操作步骤包括：
[0039]
s31：将一个二维含噪信号表示为：
[0040]
y(i,k)＝ξ(i,k)+η(i,k)
ꢀꢀꢀ
(9)，
[0041]
其中，ξ(i,k)为不含噪信号，η(i,k)为加性随机噪声，i，k分别为二维信号的行、列元素序号；
[0042]
s32：定义自适应wiener滤波所遵循的最小均方误差准则为：
[0043][0044]
其中，p、q分别为二维信号的行、列元素个数；
[0045]
s33：选取二维自适应wiener滤波的局域窗，且在局域窗中，采样点的局部均值和方差分别表示为：
[0046][0047][0048]
其中，r1×
r2为局域窗的维度；
[0049]
则对于该区域内的采样点，其二维自适应wiener滤波可表示为：
[0050][0051]
其中，ν2为噪声方差，若噪声方差未知，则ν2为所有局部估计方差的均值；
[0052]
s34：对于公式(8)所得波原子分解系数c
μ
，根据公式(11)、(12)分别求出其均值和方差，然后根据公式(13)求出二维自适应wiener滤波后的系数值
[0053]
进一步地，步骤4的具体操作步骤包括：
[0054]
s41：对滤波后的波原子系数进行重构，且重构依据的公式为：
[0055][0056]
其中，是合适的紧支系数，即为滤波后的波原子系数
[0057]
s42：根据移位不变性，将步骤1中循环平移的数据再进行反循环平移，且反循环平移所依据的公式为：
[0058][0059]
其中，是平移算子，h1，h2是网格上的平移量；
[0060]
s43：最终得到随机噪声消减、有效同相轴增强的地震记录。
[0061]
本方法与现有技术相比，具有以下有益效果：
[0062]
本发明充分利用波原子分解对具有纹理特征图像识别的优势，以及自适应wiener滤波在均方误差意义下能达到最优滤波的特点，对采集到的含噪地震记录进行循环平移及二维波原子分解，并且对分解所得系数进行自适应wiener滤波，再对滤波后的系数进行波原子重构，最后进行反循环平移得到有效信号增强的地震记录，其中背景噪声得到有力压制，同时有效同相轴得到显著增强，利于后续地震资料的解释和分析。
附图说明
[0063]
图1是含噪地震记录波原子分解系数自适应维纳滤波处理流程图；
[0064]
图2(a)-(e)是含有线性同相轴的模拟地震记录滤波实验结果图；
[0065]
图3(a)-(e)是含有双曲同相轴的模拟地震记录滤波实验结果图；
[0066]
图4(a)-(e)是实际地震记录四种方法滤波效果对比图。
具体实施方式
[0067]
为了使本领域的普通技术人员能更好地理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。
[0068]
参考附图1可以得知，基于波原子分解自适应滤波的地震勘探信号增强方法，包括以下步骤：
[0069]
步骤1：选取一幅二维含噪地震记录，对其数据进行二维循环平移；
[0070]
定义含噪地震记录为u(x1,x2)，依据公式(1)对其进行循环平移：
[0071][0072]
其中，为平移算子，h1，h2为网格上的平移量；
[0073]
对数据进行二维循环平移能够改善在图像边缘或纹理等不连续点邻域内产生伪吉布斯现象而导致有效信号失真的情况；
[0074]
步骤2：对循环平移后的含噪地震记录进行二维波原子分解，得到波原子分解系数；
[0075]
定义波原子为μ＝(j,m,n)，m＝(m1,m2)和n＝(n1,n2)，五个参量j,m1,m2,n1,n2是整数值，定义相空间中一点(x
μ
,ω
μ
)为x
μ
＝2-j
n，ω
μ
＝π2
j
m，且c1、c2是两个正的常数，位置向量x
μ
和波向量ω
μ
分别表征的空域和频域中心，则基于波包框架的波原子为：
[0076][0077][0078]
定义实值波包将波原子系数c
j,m,n
作为尺度2-j
上的抽样卷积，可得：
[0079][0080]
记hilbert变换为h，定义正交基及对偶正交基分别为：
[0081][0082][0083]
记及则形成了波原子紧支框架，且有：
[0084][0085]
根据上式可得，二维波原子分解系数为：
[0086][0087]
步骤3：对波原子分解所得系数进行自适应wiener滤波，得到滤波后的波原子系数；
[0088]
设一个二维含噪信号为：
[0089]
y(i,k)＝ξ(i,k)+η(i,k)
ꢀꢀꢀ
(9)，
[0090]
其中，ξ(i,k)为不含噪信号，η(i,k)为加性随机噪声，i，k分别为二维信号的行、列元素序号；
[0091]
并且wiener滤波的目的是得到ξ(i,k)的估值而wiener滤波所遵循的最小均方误差准则为：
[0092][0093]
其中，p、q分别为二维信号的行、列元素个数；
[0094]
当进行二维滤波时，需选取滤波掩模，也就是局域窗，在局域窗中，采样点的局部均值和方差可分别表示为：
[0095][0096][0097]
其中，r1×
r2为局域窗的维度；
[0098]
则对于该区域内的采样点，其二维自适应wiener滤波可表示为：
[0099][0100]
其中，ν2为噪声方差，若噪声方差未知，则ν2为所有局部估计方差的均值；
[0101]
求出各局部区域的估计值，最后将这些估计值合为一个整体矩阵即为最终的估计值
[0102]
对于公式(8)所得波原子分解系数c
μ
，根据公式(11)、(12)分别求出其均值和方差，然后根据公式(13)求出二维自适应wiener滤波后的系数值
[0103]
步骤4：对滤波所得的波原子系数进行重构，再进行反循环平移，最终得到有效同相轴增强的地震记录；
[0104]
首先，依据公式(14)对滤波所得的波原子系数进行重构：
[0105]
[0106]
其中，是合适的紧支系数，在此即为步骤3中求出的滤波后的波原子系数
[0107]
其次，根据移位不变性，将步骤1中循环平移的量再反循环平移回来，得到有效同相轴增强的地震记录，并且反循环平移所依据的公式为：
[0108][0109]
其中，是平移算子，h1，h2是网格上的平移量；
[0110]
最终，得到随机噪声消减、有效同相轴增强的地震记录。
[0111]
实施例
[0112]
1、仿真实验
[0113]
首先，对一幅含有线性同相轴的人工合成地震记录进行实验验证，且这幅地震记录为2ms采样，共有50道。其中包含4个直线型同相轴，每个同相轴由主频为30hz的ricker子波构成。并且附图2(a)为不含噪记录，附图2(b)是信噪比为-5db的加噪记录，附图2(c)为二维多尺度小波分解阈值滤波后重构的地震记录，附图2(d)为波原子分解系数经阈值处理后重构的地震记录，附图2(e)为波原子分解系数经二维自适应wiener滤波后重构的地震记录。
[0114]
从该记录的滤波实验中可以看出，二维多尺度小波分解阈值滤波对有效同相轴的损失较为严重，且造成同相轴边缘模糊，结果如附图2(c)所示，图中黑色椭圆框和矩形框标记出的区域为同相轴模糊与损失明显的部分；而波原子分解阈值滤波效果较前者有所改善，但记录中两边地震道的同相轴损失仍较大，结果如附图2(d)所示，黑色矩形框标记出的区域为同相轴损失明显的部分；波原子分解自适应wiener滤波效果最好，在随机噪声压制和有效信号保持方面均表现优越，结果如附图2(e)所示。图中右侧的能量棒从上至下表示能量的强弱，浅色表示信号能量最强，逐渐往下颜色越深则能量越弱。据此，可以看出附图(e)中有效信号的能量(呈亮白色)较之于附图(c)和(d)更强，即有效同相轴越呈亮白色表示能量越强，同时其黑色边缘线越清晰说明其恢复提取的效果越好，图中黑色矩形框标记出的部分，与附图(c)和(d)相应部分对比较为明显。并且经定量计算所得滤波记录的信噪比，附图2(c)-(e)，分别为7.6282db、7.8338db和8.7075db。
[0115]
对上述实验中的纯净记录加入不同强度的随机噪声得到不同信噪比的含噪记录，分别采用三种方法进行滤波实验，滤波前后地震记录的信噪比数据如表1所示(保留两位小数)，其中方法1指二维多尺度小波分解阈值滤波，方法2指二维波原子分解阈值滤波，方法3指二维波原子分解自适应wiener滤波。
[0116]
表1含噪的线性同相轴记录滤波前后信噪比
[0117][0118]
[0119]
实验中，在不同的信噪比条件下，三种方法对含噪记录的信噪比均有所提高，方法1和方法2所得记录的信噪比较为接近，方法2略高一些，而本发明提出的方法3，其所得记录信噪比最高。
[0120]
其次，对一幅含有双曲同相轴的人工合成地震记录进行滤波实验，该记录中数据的采样间隔为1ms，道间距为20m，每道ricker子波的主频为35hz，共有76道，层速度分别为2300m/s，2500m/s，2900m/s。且附图3(a)为不含噪记录，附图3(b)是信噪比为-5db的加噪记录，附图3(c)为二维多尺度小波分解阈值滤波后重构的地震记录，附图3(d)为波原子分解系数经阈值处理后重构的地震记录，附图3(e)为波原子分解系数经二维自适应wiener滤波后重构的地震记录，在附图中效果对比明显的部分采用黑色椭圆框及矩形框标出。
[0121]
附图3与附图2中有效信号的能量显示方式相同，即有效同相轴能量越强则在图中越呈亮白色，并且其黑色边缘线越清晰。通过观察三种滤波方法所得的滤波记录可知，二维多尺度小波分解阈值滤波对有效同相轴造成较大的损失，且会造成同相轴边缘模糊(如附图3(c)中椭圆框和矩形框所示出部分)；波原子分解阈值滤波对有效同相轴的损失也不容忽视，双曲同相轴拱起的部分以及记录两边的地震道损失较严重(如附图3(d)中椭圆框和矩形框所示出部分)；而附图3(e)所示的波原子分解自适应wiener滤波能够在随机噪声压制和有效信号保持方面做到很好的权衡，有效同相轴保持效果较之于前两种方法更好，图中椭圆框和矩形框标出的区域为与附图3(c)和(d)对比较为明显的部分。经过计算，滤波后，三种方法所得记录的信噪比分别为7.1616db、7.4222db和9.1073db。
[0122]
同样地，对上述实验中的纯净记录加入不同强度的随机噪声得到不同信噪比的含噪记录，分别采用三种方法进行滤波实验，滤波前后的信噪比数据如表2所示(保留两位小数)，其中方法1指二维多尺度小波分解阈值滤波，方法2指二维波原子分解阈值滤波，方法3指二维波原子分解自适应wiener滤波。
[0123]
表2含噪的双曲同相轴记录滤波前后信噪比
[0124][0125]
从表2中可以看出，在加入不同强度的随机噪声后，三种方法所得记录的信噪比均有所提高，方法1和方法2所得记录的信噪比相差不大，方法2略高一些，而本发明提出的方法3，使地震记录的信噪比提高最多。
[0126]
最后，将本发明提出的方法用于实际共炮点记录的滤波处理。该共炮点记录共有168道，1ms采样，截取记录的上半部分别采用前述的三种滤波方法以及时频峰值滤波方法进行处理。其中，附图4(a)为该实际地震记录的上半部分，图4(b)为二维多尺度小波分解阈值滤波记录，图4(c)为二维波原子分解阈值滤波记录，图4(d)为时频峰值滤波记录，图4(e)为二维波原子分解自适应wiener滤波记录。虽然时频峰值滤波方法在地震勘探信号处理方面效果较为突出，且已取得了一些值得瞩目的成果，但是该方法本身存在问题，即：滤波窗长不宜太长或太短，太长虽然能够较好地压制随机噪声，但有效信号损失较大，太短又不能
很好地压制随机噪声。对于附图4(a)所示的实际记录，应用时频峰值滤波方法时选取滤波窗长为13点，能够在随机噪声压制和有效信号保持之间达到较好的权衡。
[0127]
可以看出，附图4(b)所示的二维多尺度小波阈值滤波记录和附图4(c)所示的二维波原子分解阈值滤波记录，效果明显不理想，两者对有效同相轴造成较大程度的损失及畸变；附图4(d)所示的时频峰值滤波记录虽然背景噪声有所消减，但是有效同相轴的清晰度和连续性比较有限；而附图4(e)所示的二维波原子分解自适应wiener滤波则使有效同相轴得到显著增强，看上去更清晰、连贯。其中，附图4(a)、(d)和(e)中黑色椭圆框及矩形框标出的区域为效果对比明显的部分。
[0128]
2、结果分析
[0129]
根据以上对含有线性同相轴的人工合成地震记录，含有双曲同相轴的人工合成地震记录的实验验证结果可以得知，二维波原子分解自适应wiener滤波方法较之于二维多尺度小波分解阈值滤波方法和二维波原子分解阈值滤波方法性能更优，该方法在对随机噪声进行有力压制的同时，对有效信号的增强效果也很显著；
[0130]
另外将本发明的方法用于实际共炮点记录的滤波处理，从实验验证结果来看，二维波原子分解自适应wiener滤波方法使记录中的背景噪声大大消减，同时使有效同相轴得到显著增强，看上去更清晰、连贯；
[0131]
由于实际采集到的地震记录中含有多种噪声，其中随机噪声作为背景噪声对有效信号的影响很大，为了从背景噪声中更好地恢复提取出有效信号，采取性能优越的滤波方法十分必要。本发明提出的基于波原子分解系数自适应wiener滤波的地震勘探信号增强方法，利用波原子分解对图像纹理特征的有效识别与自适应wiener滤波在最小均方误差准则下能达到最优线性滤波的特点，使其在处理地震记录时能够很好地捕捉地震记录中的有效同相轴特征，从随机噪声环境中更好地恢复提取有效同相轴。
[0132]
本说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。尽管参照前述实施例对本发明专利进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。