CN105607125B - 基于块匹配算法和奇异值分解的地震资料噪声压制方法 - Google Patents

基于块匹配算法和奇异值分解的地震资料噪声压制方法 Download PDF

Info

Publication number
CN105607125B
CN105607125B CN201610024650.9A CN201610024650A CN105607125B CN 105607125 B CN105607125 B CN 105607125B CN 201610024650 A CN201610024650 A CN 201610024650A CN 105607125 B CN105607125 B CN 105607125B
Authority
CN
China
Prior art keywords
mrow
block
msub
sub
msubsup
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610024650.9A
Other languages
English (en)
Other versions
CN105607125A (zh
Inventor
孟可心
李娟�
李月
杜海龙
王代香
吕伟力
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN201610024650.9A priority Critical patent/CN105607125B/zh
Publication of CN105607125A publication Critical patent/CN105607125A/zh
Application granted granted Critical
Publication of CN105607125B publication Critical patent/CN105607125B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364Seismic filtering
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/30Noise handling
    • G01V2210/32Noise reduction
    • G01V2210/324Filtering

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,解决现有方法由于没有考虑到地震资料同相轴固有的局部和非局部相似性,造成的噪声压制不充分和有效信息衰减的问题,该方法包括相似性分组、三维变换奇异值滤波去噪、根据基础估计重新分组并计算维纳收缩系数和维纳滤波四个步骤。该方法在理论上采用块匹配算法,考虑到地震资料同相轴的局部和非局部相似性的固有特点,采用二维奇异值分解和一维奇异值分解结合,充分利用地震资料的结构信息,实现地震资料的稀疏表示,能够较大程度地减小有效信号衰减,并提高噪声压制效果。

Description

基于块匹配算法和奇异值分解的地震资料噪声压制方法
技术领域
本发明属于地震资料处理领域,具体涉及一种基于块匹配算法和奇异值分解的地震资料噪声压制方法。
背景技术
有关于地震资料噪声压制,并已经成熟应用到实际地震资料处理中的方法主要包括:自适应滤波,f-k滤波,多项式拟合,独立成分分析,时频峰值滤波等。M.Naghizadeh在Seismic data interpolation and denoising in the frequency-wave number domain中提出的f-k滤波,根据信号与噪声的视速度差异完成噪声压制,主要应用于去除频率较低的面波噪声。多项式拟合的拟合中心点随正交多项式系数而改变,导致处理后数据主频降低、断点变形。
随着基于多尺度分析的小波变换方法被引入到地震信号去噪中来,多分辨率分析方法为地震资料处理开辟了新的方向。Shucong L和Ergen G等在Seismic Data DenoisingSimulation Research Based on Wavelet Transform所使用的小波变换,在处理二维及更高维的信号方面存在方向局限性,并不能很好的描述有效信号的方向信息,并且二维滤波方法时常会引起信号畸变,导致出现虚假同相轴。
近年来如Ridgelet变换、Curvelet变换、Contourlet变换、Sharlet变换等具有多分辨率、多方向性、多尺度化的几何分析法成为热门方法,具体可以参考:王梅和侯振杰等的基于Ridgelet变换的多尺度去噪算法研究;孟阁阁,王德利等的基于2D Curvelet变换的多炮地震数据去噪方法研究;刘成明等的基于Shearlet变换的地震随机噪声压制。它们通过少数非零变换系数重构有效信号,在各向异性上实现了良好的稀疏表示。但上述方法采用固定基函数对信号进行分解,忽略了地震同相轴固有的时空连续性,导致噪声压制效果不佳。
发明内容
本发明的目的在于解决现有方法由于没有考虑到地震资料同相轴固有的局部和非局部相似性,造成的噪声压制不充分和有效信息衰减的问题,提出一种基于地震资料同相轴固有的局部和非局部相似性的,可以明显提高地震资料信噪比,并能较大程度保留有效信息的噪声压制方法。
本发明的目的是通过以下方案实现的:
一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,包括以下步骤:
步骤一、相似性分组:将整个地震资料切割成过完备的子块,遍历地选择参考子块,计算参考子块与一定邻域内的其他子块的欧氏距离,判断两子块是否相似,并把相似二维子块按照三维数据的方式存放起来成为相似组;
步骤二、三维奇异值滤波去噪:对相似组中每个子块先进行二维离散余弦变换,选用频域二维奇异值分解滤波,再将相似块向量化,进行频域一维奇异值分解,并进行奇异值保留去噪;二维离散余弦反变换后,将过完备子块加权聚合成原始尺寸地震资料,得到基础估计;
步骤三、计算维纳收缩系数:将基础估计重新分块分组,并且通过新的相似组计算变换域的维纳收缩系数;
步骤四、维纳滤波:按照新的相似分组,对原始信号进行三维离散余弦变换,采用维纳滤波实现噪声压制,三维离散余弦反变换并聚合重构得到最终噪声压制后的地震资料。
所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其中,步骤一相似性分组包括以下具体步骤:
1)预处理:先将整个地震资料按照一定的步进长度过完备的分割成N1×N1大小的时-空域子块Yx;然后对各个子块的二维频谱进行阈值r′处理,将大于阈值的频谱值保留下来,小于阈值的频谱值置0;
2)计算块间距离:将预处理后的子块频谱r′(T(Yx))作为计算块间距离的元素,选定一个参考块YR,选定搜索窗大小为Nh×Nh,通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离:
3)堆叠成三维数组:当两个子块的距离d(YR,Yx)小于阈值τ时,判定两块相似,然后将所有与同一参考块相似的二维子块,连同参考块一起存放成一个三维数组,即相似组ZR
所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其中,步骤二三维奇异值滤波去噪包括以下具体步骤:
1)二维离散余弦变换:进行二维离散余弦变换,将相似组中各个子块变换到频域;
2)奇异值处理:
(1)二维奇异值分解:对于一个相似组中的任一二维子块A,A是N1×N1的,它的特征矩阵U2D和V2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得,将Ai,i∈[0,N1-1]投影到U2D和V2D上,得到投影系数Λi,满足:
Λi即奇异值矩阵,将奇异值矩阵中大于阈值T1的奇异值保留,小于T1的奇异值置0,然后重构;
(2)一维奇异值分解:将相似组中所有子块向量化,进行一维奇异值分解,获得奇异值矩阵并进行重构;
3)离散余弦反变换、加权聚合:经所述步骤2)处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块,将其进行二维离散余弦反变换,然后对二维离散余弦反变换结果进行加权重构,得到实现大部分噪声压制的该相似组所有子块对应位置的时域的的基础估计
其中,NR是未置0的系数个数
所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其中,步骤三计算维纳收缩系数具体包括以下步骤:
1)重新分块:将经过所述步骤二获得的基础估计过完备地分成N1×N1的子块选择参考块重新计算子块间欧氏距离:
其中,是经过所述步骤二得到的基础估计;
2)堆叠成三维数组:当两个子块的距离小于阈值τ时,判定两块相似,并将分组方法传递给原始含噪地震资料,将判定为相似的子块作为一个数组存放起来形成一个新的相似组
3)三维离散余弦变换:首先对相似组中每个二维块独立地进行二维离散余弦变换,然后对所有子块的同一位置的数据点进行一维离散余弦变换,三维离散余弦变换后,地震资料被变换到频域;
4)计算维纳收缩系数:
根据基础估计计算维纳收缩系数WR
其中,σ2是噪声方差,T3D为三维离散余弦变换,的分组块。
所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其中,步骤四维纳滤波具体包括以下步骤:
1)维纳滤波:在所述步骤三中新的相似组已经经过三维离散余弦变换,在此将频域块与对应的维纳收缩系数相乘实现滤波;维纳滤波的计算方法是:
其中,σ2是噪声方差,T3D为三维离散余弦变换,是经过所述步骤二获得的基础估计的分组块;
2)三维离散余弦反变换:按照先一维离散余弦反变换,再二维离散余弦反变换的倒序方式,对每个子块依次实现反变换;
3)加权聚合:对二维离散余弦反变换结果进行加权聚合,得到最终噪声压制后的地震资料
其中,
与现有技术方案相比,本发明的方法在理论上采用块匹配算法,考虑到地震资料的固有特点——同相轴的局部和非局部相似性;采用二维奇异值分解和一维奇异值分解结合,充分利用地震资料的结构信息,实现地震资料的稀疏表示,能够较大程度地减小有效信号衰减,并提高噪声压制效果。从数据方面来看,通过仿真模型可以看出,表征信号保留程度的均方根误差减小了0.002,表征噪声压制效果的信噪比增加了大约3.4db,在250~500Hz的高频段噪声压制效果尤其明显。
附图说明
图1为本发明方法流程图。
图2为块匹配图示。
图3为纯净地震仿真信号
图4为图3所示仿真信号的加入-5db噪声地震仿真信号
图5为原始方法处理后的地震仿真信号
图6为本发明方法处理后的地震仿真信号
图7为第30道仿真信号处理对比图
图8为第30道仿真信号放大对比图
图9为第30道信号幅度谱对比图
图10为实际带噪地震资料
图11为原始方法处理后的实际地震资料
图12为本发明方法处理后的实际地震资料
图13为实际信号放大对比图
具体实施方式
本发明提出的地震资料噪声压制方法,是一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,该算法结合了三维块匹配算法和奇异值分解算法。
本发明中块匹配是指将待处理信号过完备的分成特定大小的子块,遍历性地在一定搜索邻域内寻找满足某种相似规范的相似块,然后将所有与同一参考块匹配的二维块与参考块一起,按照三维数组的形式存放起来。过完备是指分成的子块是相互重叠的,其中心点数据的序号差为分块时选择的步进长度,如图2所示。遍历性地搜索是指按顺序将每一个子块作为参考块,考虑后续的计算过程。相似性分组的处理方法充分的考虑到同相轴的非局部相似性。
在频率域,不同道的地震信号具有相同的有效频带,并且在地震道的方向上这些频率具有可预测性、横相干性;而噪声是随机的,频带较宽,也不存在可预测性和相关性。信号在频率域能得到比时域更稀疏的表示方法,由此,本发明采用在频率域进行奇异值分解去噪。
考虑到从属于不同块的同相轴间具有非局部相似性,本发明采用一维奇异值分解实现同相轴非局部冗余的去除。一维奇异值分解充分保留了地震资料的非局部结构信息,在实现噪声压制的同时,也能较大程度地保留有效信息。
考虑到从属于同一块的同相轴间具有局部相似性,本发明采用二维奇异值分解实现同相轴局部冗余的去除。二维奇异值分解充分保留了地震资料的局部结构信息,在实现噪声压制的同时,也较大程度地保留了有效信息。
三维块匹配算法和奇异值分解相结合充分考虑了同相轴局部和非局部的相似性,避免了现有方法中没有考虑同相轴局部和非局部的相似性而导致的噪声压制不充分和信号畸变。
本发明提出的基于块匹配算法和奇异值分解的地震资料噪声压制方法主要分为4步,包括相似性分组、三维变换奇异值滤波去噪、根据基础估计重新分组并计算维纳收缩系数和维纳滤波四个步骤。如图1所示。
步骤一、相似性分组:将整个地震资料切割成过完备的子块,遍历地选择参考子块,计算参考子块与一定邻域内的其他子块的欧氏距离,判断两子块是否相似,并把相似二维子块按照三维数据的方式存放起来成为三维数组。
步骤二、三维奇异值滤波去噪:对相似组中每个子块先进行二维离散余弦变换,选用频域二维奇异值分解滤波。再将相似块向量化,进行频域一维奇异值分解,进行奇异值保留去噪。二维离散余弦反变换后,将过完备子块加权聚合成原始尺寸地震资料,得到基础估计。
步骤三、计算维纳收缩系数:将基础估计重新分块分组,并且通过新的相似分组计算变换域的维纳收缩系数。
步骤四、维纳滤波:按照新的相似分组,对原始信号进行三维离散余弦变换,采用维纳滤波实现噪声压制,三维离散余弦反变换并聚合重构得到最终估计,即最终噪声压制后的地震资料。
下面对各步骤进一步介绍如下:
步骤一、相似性分组:
1.预处理
先将整个地震资料按照一定的步进长度过完备的分割成N1×N1大小的时-空域子块Yx,准备计算其相似性差值。但是,当地震资料中包含噪声较多时,可能由于噪声方差过大而引起相似规范的计算误差,从而引起分组错误。因此,在计算块间距离之前,首先对各个子块的二维频谱进行阈值r′处理。将大于阈值的频谱值保留下来,小于阈值的频谱值置0。预处理过程能够减小下一步计算块间距离过程中噪声对相似判断的错误影响。
2.计算块间距离
相似性的规范一般采用计算信号两个块之间的一些距离参数,例如欧几里得距离。本发明将预处理后的子块频谱r′(T(Yx))作为计算块间距离的元素,选定一个参考块YR,选定搜索窗大小为Nh×Nh,通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离:
通过计算块间距离能够得知参考块邻域中的其他子块与这一参考块的相似程度,为下一步相似性分组做好了准备。
3.堆叠成三维数组
通过计算块间距离,已经得知参考块的一定邻域内其他块与当前参考块的相似程度,为了后续滤波过程的顺利实施,需要将判定为相似的时-空域子块作为一个数组存放起来。当两个子块的距离d(YR,Yx)小于阈值τ时,判定两块相似,然后将所有与同一参考块相似的二维子块,连同参考块一起存放成一个三维数组,即相似组ZR。经过分组后,具有非局部相似性的子块被存放在一起,为一维奇异值分解做好了准备。
步骤二、三维奇异值滤波去噪
1.二维离散余弦变换
对相似组ZR进行三维奇异值滤波之前,首先进行二维离散余弦变换,将相似组ZR中各个子块变换到频域。上面提到过在频率域,不同道的地震信号具有相同的有效频带,并且在地震道的方向上这些频率具有可预测性、横相干性;而噪声是随机的,频带较宽,也不存在可预测性和相关性。也就是说在频域中,有效信号的能量能更集中在某几个少数不为0的点上,即稀疏表示。这样做扩大了有效信息与噪声的差异,能够增强奇异值滤波的滤波效果。
2.奇异值处理
(1)二维奇异值分解
经过二维离散余弦变换,子块的地震数据变换到频域,实现初步的稀疏表示。对于一个相似组中的任一二维子块A,A是N1×N1的,它的特征矩阵U2D和V2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得,将Ai,i∈[0,N1-1]投影到U2D和V2D上,得到投影系数Λi,满足
Λi即奇异值矩阵。奇异值分解使得地震资料能用更少的非零系数表示完全,实现了进一步稀疏表示。而且,奇异值矩阵中较大的奇异值代表地震资料中的有效信息,较小的奇异值代表噪声信息。因此,将大于阈值T1的奇异值保留,小于T1的奇异值置0,然后重构,这一步考虑了局部相似性,实现了块内冗余的去除。
(2)一维奇异值分解
二维奇异值分解只考虑到地震资料的局部相似性,而没有考虑到其非局部相似性。因此,将相似组中所有子块向量化,即将二维数据子块变成一个列向量,进行一维奇异值分解。
每个子块是N1×N1的,假设此参考块有M-1个相似块,那么,构造出来的相似矩阵B的尺寸为N1×N1行,M列。对于此矩阵必然存在两个正交矩阵满足
Λ1D即奇异值矩阵,同样将Λ1D中将大于阈值T2的奇异值保留,小于T2的奇异值置0,然后重构,这一步骤考虑了同相轴的非局部相似性,实现了块间冗余的去除。
3.离散余弦反变换、加权聚合
处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块,将其进行二维离散余弦反变换,得到该相似组所有子块对应位置的时域的基础估计但块是相互重叠的,因此需要加权重构。对应权重计算方法是:
由于噪声方差σ2和未置0的系数个数NR的存在,权重代表了该参考块中有效信息的占比,有效信息越多,权重越大。相互重叠的过完备块的基础估计按照下式方法聚合。
由此,得到已经实现大部分噪声压制的基础估计
步骤三、计算维纳收缩系数
1.重新分块
由于基础估计已经实现大部分的噪声压制,能够实现更精确的相似性分组。将基础估计过完备地分成N1×N1的子块选择参考块按照同样的方法重新计算子块间欧氏距离:
这里是上一步得到的基础估计,得到的是基础估计中每一参考块与其周围邻域中各个子块的相似性差值,为下一步相似性分组做好准备。
2.堆叠成三维数组
通过计算块间距离,已经得知参考块的一定邻域内其他块与当前参考块的相似程度,为了后续滤波过程的顺利实施,需要将判定为相似的子块作为一个数组存放起来。当两个子块的距离小于阈值τ时,判定两块相似,并将分组方法传递给原始含噪地震资料,形成一个新的三维数组重新分组后得到更为精确的非局部相似性,能更为准确地保留地震资料的非局部结构信息。
3.三维离散余弦变换
在进行滤波处理之前,首先将地震资料由时域变换到频域,本发明采用的是三维离散余弦变换。三维离散余弦变换是首先对相似组中每个二维块独立地进行二维离散余弦变换,然后对所有块的同一位置的数据点进行一维离散余弦变换。三维离散余弦变换后,地震资料被变换到频域,实现了初步的稀疏表示。
4.计算维纳收缩系数
另外,根据基础估计计算下一步滤波所需要的维纳收缩系数WR
这里σ2是噪声方差,T3D为三维离散余弦变换,是基础估计的分组块。通过维纳系数的计算方式可以看出,维纳系数代表了该子块的频域能量与估计的含噪信息的频域能量的比值,即有效信息越多,维纳收缩系数越大,这也是它的去噪原理。
步骤四、维纳滤波
1.维纳滤波
在第三个大步骤中,通过各个小步的计算,已经完成了所有维纳滤波的准备过程,在这一部分中将实现维纳滤波去除地震资料中的噪声。维纳滤波的计算方法是:
在上述步骤三中分组块已经经过三维离散余弦变换,在此需要将频域块与对应的维纳收缩系数相乘实现滤波。上面已经提到包含有效信息的数据点对应的维纳收缩系数较大,也就是说通过维纳滤波子块中代表有效信息的数据点得以保留,而代表噪声的数据点被压制,实现了地震资料子块的噪声去除。
2.三维离散余弦反变换
进行三维离散余弦反变换将频域子块变回到时域,得到过完备的块的估计。三维离散余弦反变换过程按照先一维离散余弦反变换,再二维离散余弦反变换的倒序方式,对每个子块依次实现反变换。但此时块的估计是相互重叠的,需要下一步的加权聚合得到最终的噪声压制结果。
3.加权聚合
同样,对过完备块进行加权求和,权重是根据维纳收缩系数计算得来的。
权重同样代表该块中有效信息的占比。聚合方法与之前类似:
这样就得到了最终噪声压制后的地震资料
实施例
一.仿真信号模型实施例
仿真条件如下:在Matlab实验平台,采用100道,每道2000个采样点的合成地震信号仿真模型来进行算法的性能验证。该记录包含三条同相轴,主频分别为45Hz、35Hz和25Hz,视速度分别为1000m/s、1400m/s和1500m/s,采样频率是1000Hz,炮检距为10m。不包含噪声的纯净地震信号仿真模型如图3所示。
为验证本发明方法对低信噪比地震资料的噪声压制效果,在纯净地震信号仿真模型中加入-5db的高斯白噪声,如图4所示,可以看到有效信息大量淹没在噪声之中。
1、相似性分组
将2000×100的带噪地震仿真模型按照块步进长度为2,分成1000×50个8×8的子块。首先将每个子块通过二维傅里叶变换转换到频域,经过阈值为0.9的预处理。频域值大于0.9的数据点保留,小于0.9的数据点置零。然后以每个子块为中心,在其周围19×19的邻域内搜索其相似块,相似规范为2-范数距离。把2-范数距离小于3的判为两块相似,与其对应的参考块存放在一起作为一个相似组,而距离大于3的则不做处理。
2、三维奇异值滤波去噪
将时域同一相似组中的每个子块进行二维离散余弦变换转换到频域,实现初步的稀疏表示。然后对每个子块独立的进行二维奇异值分解,以1.2作为阈值处理二维奇异值,大于1.2的奇异值保留下来,小于1.2的奇异值作为噪声分量置0,然后重构回去,得到该相似组中每个子块的估计,是去除块内冗余的结果。然后将所有子块向量化,即将二维数据块转化为一个列向量,然后依次按列存放成一个二维数据块,做一维奇异值分解。这一步以1作为阈值,重构得到的是去除这一相似组中块间冗余的结果。由于子块是相互重叠的,按照上文提到的公式计算权重,加权求和即得到整个地震仿真模型的基础估计。
3、计算维纳收缩系数
将2000×100的地震资料基础估计按照块步进长度为2,分成1000×50个8×8的子块。然后以每个子块为中心,在其周围19×19的邻域内搜索其相似块,相似规范为2-范数距离。把2-范数距离小于1的判为两块相似,与其对应的参考块存放在一起作为一个相似组,而距离大于1的则不做处理。同时将这种分组方法,通过子块序号的方式传递给原始带噪地震资料,将原始带噪地震资料重新分组。
将原始带噪地震资料新的相似组中的每个二维块进行二维离散余弦变换,然后再对每个二维块相同位置的数据点进行一维离散余弦变换,将时域信号变换到频域,实现地震资料的初步稀疏表示。并通过各数据块的频谱计算对应的维纳收缩系数,为下一步滤波做好准备工作。
4、维纳滤波
将每一子块与其对应的维纳收缩系数相乘实现滤波。这一步不需要阈值处理,因为维纳收缩系数的计算方法使得包含有效信息的点系数较大,噪声点系数较小。因此子块与其相乘就能实现有效信息的保留和噪声压制效果。滤波后,将频域子块按照先一维离散余弦反变换,再二维离散余弦反变换的顺序变换回时域。然后将重叠的子块估计加权聚合,就能得到2000×100的地震资料仿真模型最终的噪声压制效果。
经过原始三维块匹配算法处理后的噪声压制效果如图5所示,可以看到与图3的纯净信号相比,图5仍包含少量噪声。而图6,本发明算法处理后的地震信号明显剩余噪声更少,有更强的噪声压制效果。表1从均方根误差和信噪比两个方面,量化地证明本发明算法取得的更好的效果。
表1去噪效果对比表
均方根误差 信噪比
原始三维块匹配 0.025 13.43db
本发明算法 0.023 16.86db
图7在100道地震信号仿真模型中选取了第30道地震信号,图8放大其中一条同相轴进行观察。可以清晰地看到本发明算法处理后的信号曲线在不含信号部分更为平滑,包含信号的部分曲线与纯净信号曲线更为贴合。这就证明本发明算法有更好的噪声压制效果,同时也能更好地保留有效信号的结构信息。
图9从第30道地震信号的幅度谱的角度证明这一结论。从幅度谱来看,本发明算法在250Hz~500Hz的频率段更为贴合纯净信号的幅度谱,而原始三维块匹配在这一频段仍包含噪声。
二.实际地震资料实施例
实际地震资料是一张168道,每道6000点的共炮点记录,如图10。采样频率为1000Hz,炮检距为30m。该资料中包含随机噪声、较为混乱的同相轴和设备带来的其他噪声。
1、相似性分组
将6000×168的带噪地震仿真模型按照块步进长度为4,分成1500×42个8×8的子块。首先将每个子块通过二维傅里叶变换转换到频域,经过阈值为1×10-5的预处理。频域值大于1×10-5的数据点保留,小于1×10-5的数据点置零。然后以每个子块为中心,在其周围39×39的邻域内搜索其相似块,相似规范为2-范数距离。把2-范数距离小于2×10-9的判为两块相似,与其对应的参考块存放在一起作为一个相似组,而距离大于2×10-9的则不做处理。
2、三维奇异值滤波去噪
将时域同一相似组中的每个子块进行二维离散余弦变换转换到频域,实现初步的稀疏表示。然后对每个子块独立的进行二维奇异值分解,以1×10-4作为阈值处理二维奇异值,大于1×10-4的奇异值保留下来,小于1×10-4的奇异值作为噪声分量置0,然后重构回去,得到该相似组中每个子块的估计,是去除块内冗余的结果。然后将所有子块向量化,即将二维数据块转化为一个列向量,然后依次按列存放成一个二维数据块,做一维奇异值分解。这一步以0.0012作为阈值,重构得到的是去除这一相似组中块间冗余的结果。由于子块是相互重叠的,按照上文提到的公式计算权重,加权求和即得到整个地震仿真模型的基础估计。
3、计算维纳收缩系数
将6000×168的地震资料基础估计按照块步进长度为4,分成1500×42个8×8的子块。然后以每个子块为中心,在其周围39×39的邻域内搜索其相似块,相似规范为2-范数距离。把2-范数距离小于1×10-15的判为两块相似,与其对应的参考块存放在一起作为一个相似组,而距离大于1×10-15的则不做处理。同时将这种分组方法,通过子块序号的方式传递给原始带噪地震资料,将原始带噪地震资料重新分组。
将原始带噪地震资料新的相似组中的每个二维块进行二维离散余弦变换,然后再对每个二维块相同位置的数据点进行一维离散余弦变换,将时域信号变换到频域,实现地震资料的初步稀疏表示。并通过各数据块的频谱计算对应的维纳收缩系数,为下一步滤波做好准备工作。
4、维纳滤波
将每一子块与其对应的维纳收缩系数相乘实现滤波。这一步不需要阈值处理,因为维纳收缩系数的计算方法使得包含有效信息的点系数较大,噪声点系数较小。因此子块与其相乘就能实现有效信息的保留和噪声压制效果。滤波后,将频域子块按照先一维离散余弦反变换,再二维离散余弦反变换的顺序变换回时域。然后将重叠的子块估计加权聚合,就能得到6000×168的地震资料仿真模型最终的噪声压制效果。
原始三维块匹配和本发明算法处理后的效果分别如图11和图12所示。可以清晰看到本发明算法处理后的地震资料有更为清晰的同相轴。放大后观察,也能看出图13(b)本算法处理后含的噪声比图13(a)原始三维块匹配算法处理后的噪声更少。

Claims (5)

1.一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其特征在于,包括以下步骤:
步骤一、相似性分组:将整个地震资料切割成过完备的子块,遍历地选择参考子块,计算参考子块与一定邻域内的其他子块的欧氏距离,判断两子块是否相似,并把相似二维子块按照三维数据的方式存放起来成为相似组;
步骤二、三维奇异值滤波去噪:对相似组中每个子块先进行二维离散余弦变换,选用频域二维奇异值分解滤波,再将相似块向量化,进行频域一维奇异值分解,并进行奇异值保留去噪;二维离散余弦反变换后,将过完备子块加权聚合成原始尺寸地震资料,得到基础估计;
步骤三、计算维纳收缩系数:将基础估计重新分块分组,并且通过新的相似组计算变换域的维纳收缩系数;
步骤四、维纳滤波:按照新的相似分组,对原始信号进行三维离散余弦变换,采用维纳滤波实现噪声压制,三维离散余弦反变换并聚合重构得到最终噪声压制后的地震资料。
2.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其特征在于,所述步骤一相似性分组包括以下具体步骤:
1)预处理:先将整个地震资料按照一定的步进长度过完备的分割成N1×N1大小的时-空域子块Yx;然后对各个子块的二维频谱进行阈值r′处理,将大于阈值的频谱值保留下来,小于阈值的频谱值置0;
2)计算块间距离:将预处理后的子块频谱r′(T(Yx))作为计算块间距离的元素,选定一个参考块YR,选定搜索窗大小为Nh×Nh,通过2-范数的计算参考块与其相应搜索窗内其他子块的块间距离:
<mrow> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> <mo>,</mo> <msub> <mi>Y</mi> <mi>x</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>T</mi> <mo>(</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> <mo>-</mo> <msup> <mi>r</mi> <mo>&amp;prime;</mo> </msup> <mrow> <mo>(</mo> <mi>T</mi> <mo>(</mo> <msub> <mi>Y</mi> <mi>x</mi> </msub> <mo>)</mo> </mrow> <mo>)</mo> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> </mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>N</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> </mrow>
3)堆叠成三维数组:当两个子块的距离d(YR,Yx)小于阈值τ时,判定两块相似,然后将所有与同一参考块相似的二维子块,连同参考块一起存放成一个三维数组,即相似组ZR
3.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其特征在于,所述步骤二三维奇异值滤波去噪包括以下具体步骤:
1)二维离散余弦变换:进行二维离散余弦变换,将相似组中各个子块变换到频域;
2)奇异值处理:
(1)二维奇异值分解:对于一个相似组中的任一二维子块A,A是N1×N1的,它的特征矩阵U2D和V2D分别根据它的行-行协方差矩阵和列-列协方差矩阵求得,将Ai,i∈[0,N1-1]投影到U2D和V2D上,得到投影系数Λi,满足:
<mrow> <msub> <mi>A</mi> <mi>i</mi> </msub> <mo>=</mo> <msub> <mi>U</mi> <mrow> <mn>2</mn> <mi>D</mi> </mrow> </msub> <msub> <mi>&amp;Lambda;</mi> <mi>i</mi> </msub> <msubsup> <mi>V</mi> <mrow> <mn>2</mn> <mi>D</mi> </mrow> <mi>T</mi> </msubsup> </mrow>
其中,Ai为所述二维子块矩阵A的第i列;Λi即奇异值矩阵,将奇异值矩阵中大于阈值T1的奇异值保留,小于T1的奇异值置0,然后重构;
(2)一维奇异值分解:将相似组中所有子块向量化,进行一维奇异值分解,获得奇异值矩阵并进行重构;
3)离散余弦反变换、加权聚合:经所述步骤2)处理后的奇异值矩阵重构得到的是噪声压制后的频域数据块,将其进行二维离散余弦反变换,然后对二维离散余弦反变换结果进行加权重构,得到实现大部分噪声压制的该相似组所有子块对应位置的时域的基础估计
<mrow> <msup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msup> <mo>=</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> </munder> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Z</mi> <mi>R</mi> </msub> </munder> <msubsup> <mi>&amp;omega;</mi> <mi>R</mi> <mrow> <mi>h</mi> <mi>t</mi> </mrow> </msubsup> <mo>&amp;times;</mo> <msub> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> </msub> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> </munder> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Z</mi> <mi>R</mi> </msub> </munder> <msubsup> <mi>&amp;omega;</mi> <mi>R</mi> <mrow> <mi>h</mi> <mi>t</mi> </mrow> </msubsup> </mrow> </mfrac> </mrow>
其中,NR是未置0的系数个数。
4.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其特征在于,所述步骤三计算维纳收缩系数具体包括以下步骤:
1)重新分块:将经过所述步骤二获得的基础估计过完备地分成N1×N1的子块选择参考块重新计算子块间欧氏距离:
<mrow> <mi>d</mi> <mrow> <mo>(</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>,</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>x</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <mo>|</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>-</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>x</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> </mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>N</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mfrac> </mrow>
其中,是经过所述步骤二得到的基础估计;
2)堆叠成三维数组:当两个子块的距离小于阈值τ时,判定两块相似,并将分组方法传递给原始含噪地震资料,将判定为相似的子块作为一个数组存放起来形成一个新的相似组
3)三维离散余弦变换:首先对相似组中每个二维块独立地进行二维离散余弦变换,然后对所有子块的同一位置的数据点进行一维离散余弦变换,三维离散余弦变换后,地震资料被变换到频域;
4)计算维纳收缩系数:
根据基础估计计算维纳收缩系数WR
<mrow> <msub> <mi>W</mi> <mi>R</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <msub> <mi>T</mi> <mrow> <mn>3</mn> <mi>D</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>|</mo> </mrow> <mrow> <mo>|</mo> <msub> <mi>T</mi> <mrow> <mn>3</mn> <mi>D</mi> </mrow> </msub> <mrow> <mo>(</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> <mrow> <mi>b</mi> <mi>a</mi> <mi>s</mi> <mi>i</mi> <mi>c</mi> </mrow> </msubsup> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <msup> <mi>&amp;sigma;</mi> <mn>2</mn> </msup> </mrow> </mfrac> </mrow>
其中,σ2是噪声方差,T3D为三维离散余弦变换,的分组块。
5.如权利要求1所述的一种基于块匹配算法和奇异值分解的地震资料噪声压制方法,其特征在于,所述步骤四维纳滤波具体包括以下步骤:
1)维纳滤波:在所述步骤三中新的相似组已经经过三维离散余弦变换,在此将频域块与对应的维纳收缩系数相乘实现滤波;维纳滤波的计算方法是:
其中,σ2是噪声方差,T3D为三维离散余弦变换,是经过所述步骤二获得的基础估计的分组块;
2)三维离散余弦反变换:按照先一维离散余弦反变换,再二维离散余弦反变换的倒序方式,对每个子块依次实现反变换;
3)加权聚合:对二维离散余弦反变换结果进行加权聚合,得到最终噪声压制后的地震资料
<mrow> <msup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msup> <mo>=</mo> <mfrac> <mrow> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> </munder> <munder> <mo>&amp;Sigma;</mo> <msubsup> <mi>Z</mi> <mi>R</mi> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msubsup> </munder> <msubsup> <mi>&amp;omega;</mi> <mi>R</mi> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msubsup> <mo>&amp;times;</mo> <msubsup> <mover> <mi>Y</mi> <mo>^</mo> </mover> <mi>R</mi> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msubsup> </mrow> <mrow> <munder> <mo>&amp;Sigma;</mo> <msub> <mi>Y</mi> <mi>R</mi> </msub> </munder> <munder> <mo>&amp;Sigma;</mo> <msubsup> <mi>Z</mi> <mi>R</mi> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msubsup> </munder> <msubsup> <mi>&amp;omega;</mi> <mi>R</mi> <mrow> <mi>w</mi> <mi>i</mi> <mi>e</mi> </mrow> </msubsup> </mrow> </mfrac> </mrow>
其中,
CN201610024650.9A 2016-01-15 2016-01-15 基于块匹配算法和奇异值分解的地震资料噪声压制方法 Expired - Fee Related CN105607125B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610024650.9A CN105607125B (zh) 2016-01-15 2016-01-15 基于块匹配算法和奇异值分解的地震资料噪声压制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610024650.9A CN105607125B (zh) 2016-01-15 2016-01-15 基于块匹配算法和奇异值分解的地震资料噪声压制方法

Publications (2)

Publication Number Publication Date
CN105607125A CN105607125A (zh) 2016-05-25
CN105607125B true CN105607125B (zh) 2018-04-13

Family

ID=55987201

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610024650.9A Expired - Fee Related CN105607125B (zh) 2016-01-15 2016-01-15 基于块匹配算法和奇异值分解的地震资料噪声压制方法

Country Status (1)

Country Link
CN (1) CN105607125B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108037533B (zh) * 2018-01-16 2019-05-31 吉林大学 一种基于块重组期望对数似然的地震勘探噪声压制方法
CN108645920B (zh) * 2018-04-09 2020-12-22 华南理工大学 一种基于去噪和对齐的钢轨超声探伤的直达波抑制方法
CN108680950B (zh) * 2018-05-16 2019-07-26 吉林大学 一种基于自适应块匹配的沙漠地震信号位置检测方法
CN108957550B (zh) * 2018-06-28 2020-01-03 吉林大学 基于svd-ica的tsp强工业电干扰压制方法
CN109493929B (zh) * 2018-09-20 2022-03-15 北京工业大学 基于分组变量的低冗余特征选择方法
CN109171706B (zh) * 2018-09-30 2020-12-29 南京信息工程大学 基于分类匹配和分数阶扩散的心电信号去噪方法和系统
CN109991657B (zh) * 2018-11-15 2021-10-15 成都理工大学 基于逆二分递推奇异值分解的地震资料高分辨率处理方法
CN111324598A (zh) * 2018-12-17 2020-06-23 中国石油天然气股份有限公司 微震记录的去噪方法及装置
CN109633751A (zh) * 2018-12-26 2019-04-16 北京化工大学 一种基于多分辨奇异值分解的地震数据降噪方法
CN110068865B (zh) * 2019-05-09 2021-02-23 吉林大学 一种低秩矩阵逼近的沙漠地震噪声压制方法
CN112649863A (zh) * 2019-10-12 2021-04-13 中国石油化工股份有限公司 分频地震属性数据优化方法及系统
CN113138409A (zh) * 2020-01-19 2021-07-20 中国石油天然气集团有限公司 一种三维叠后地震数据处理方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101852866A (zh) * 2010-05-13 2010-10-06 中国石油天然气股份有限公司 一种叠后地震数据滤波方法
WO2013074720A1 (en) * 2011-11-18 2013-05-23 Geco Technology B.V. Noise removal from 3d seismic representation
GB2500302A (en) * 2012-01-31 2013-09-18 Cggveritas Services Sa Removing non-orthogonal coherent noise from seismic data
CN104159003A (zh) * 2014-08-21 2014-11-19 武汉大学 一种基于3d协同滤波与低秩矩阵重建的视频去噪方法及系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101852866A (zh) * 2010-05-13 2010-10-06 中国石油天然气股份有限公司 一种叠后地震数据滤波方法
WO2013074720A1 (en) * 2011-11-18 2013-05-23 Geco Technology B.V. Noise removal from 3d seismic representation
GB2500302A (en) * 2012-01-31 2013-09-18 Cggveritas Services Sa Removing non-orthogonal coherent noise from seismic data
CN104159003A (zh) * 2014-08-21 2014-11-19 武汉大学 一种基于3d协同滤波与低秩矩阵重建的视频去噪方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于频域正则维纳滤波的地震随机噪声压制;田雅男 等;《吉林大学学报(工学版)》;20151130;第45卷(第6期);第2043-2048页 *
频域奇异值分解(SVD)地震波场去噪;沈鸿雁 等;《石油地球物理勘探》;20100430;第45卷(第2期);第185-189页 *

Also Published As

Publication number Publication date
CN105607125A (zh) 2016-05-25

Similar Documents

Publication Publication Date Title
CN105607125B (zh) 基于块匹配算法和奇异值分解的地震资料噪声压制方法
CN102854533B (zh) 一种基于波场分离原理提高地震资料信噪比的去噪方法
Kingsbury et al. Wavelet transforms in image processing
CN101303764B (zh) 基于非下采样轮廓波的多传感器图像自适应融合方法
CN102158637B (zh) 基于Surfacelet变换域的空间自适应阈值视频去噪方法
Zhang et al. A patch based denoising method using deep convolutional neural network for seismic image
CN107192553A (zh) 基于盲源分离的齿轮箱复合故障诊断方法
CN105913393A (zh) 一种自适应小波阈值图像去噪算法及装置
CN105974468B (zh) 一种能够同时进行五维地震数据重建和噪声压制的方法
CN107144879B (zh) 一种基于自适应滤波与小波变换结合的地震波降噪方法
CN108805840A (zh) 图像去噪的方法、装置、终端及计算机可读存储介质
CN106771905B (zh) 一种适用于高频电流局部放电检测的脉冲提取方法
CN109164483A (zh) 多分量地震数据矢量去噪方法及多分量地震数据矢量去噪装置
CN101729157A (zh) 一种强噪声环境下的振动信号盲源分离方法
CN104780127B (zh) 一种时延‑多普勒r‑l解卷多路径信道估计方法
CN103208097A (zh) 图像多方向形态结构分组的主分量分析协同滤波方法
CN107251053A (zh) 一种降低有损压缩图像的压缩失真的方法及装置
CN104280776B (zh) 一种自适应小波阈值求取方法
CN113204051B (zh) 一种基于变分模态分解的低秩张量地震数据去噪方法
CN117111000A (zh) 一种基于双通道注意力残差网络的sar梳状谱干扰抑制方法
CN103784164B (zh) 超声信号的预处理方法及系统
CN103645504A (zh) 基于广义瞬时相位及p范数负模的地震弱信号处理方法
CN106950597B (zh) 基于三边滤波的混合震源数据分离方法
CN102509268B (zh) 基于免疫克隆选择的非下采样轮廓波域图像去噪方法
CN109917459A (zh) 一种压制地震噪声的方法、装置及系统

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20180413

Termination date: 20190115

CF01 Termination of patent right due to non-payment of annual fee