基于自适应K-SVD的能量泄漏恢复研究

doi:10.11720/wtyht.2020.1229

基于自适应K-SVD的能量泄漏恢复研究

何真¹, 曹思远²^,³, 郝婳婕⁴, 段绪远⁵, 叶力夏提⁵, 张韵佩¹

1. 中国石油新疆油田分公司勘探开发研究院,新疆克拉玛依 834000

2. 中国石油大学(北京) 地球物理学院,北京 102249

3. 中国石油大学(北京) CNPC物探重点实验室,北京 102249

4. 中国石油新疆油田分公司数据公司信息研究所,新疆克拉玛依 834000

5. 中国石油新疆油田分公司新港公司,新疆克拉玛依 834000

Research on energy leakage recovery of adaptive K-SVD

HE Zhen¹, CAO Si-Yuan²^,³, HAO Hua-Jie⁴, DUAN Xu-Yuan⁵, ⁵, ZHANG Yun-Pei¹

1. Exploration and Development Research Institute of PetroChina Xinjiang Oilfield Company,Karamay 834000,China

2. Institute of Geophysics,China University of Petroleum (Beijing) ,Beijing 102249,China

3. CNPC Geophysical Key Laboratory,China University of Petroleum (Beijing),Beijing 102249,China

4. Data Company Information Institute of PetroChina Xinjiang Oilfield Company,Karamay 834000,China

5. Xingang Branch Company of PetroChina Xinjiang Oilfield Company,Karamay 834000,China

责任编辑: 叶佩

收稿日期: 2019-06-21 修回日期: 2020-01-6 网络出版日期: 2020-04-20

基金资助:

国家自然科学基金项目.  41674128
中国石油科技创新基金项目“陆上层间多次波压制技术”.  2017D-5007-0302
国家重点研发计划项目.  SQ2017YFGX030021

Received: 2019-06-21 Revised: 2020-01-6 Online: 2020-04-20

作者简介 About authors

何真(1992-),男,硕士研究生学历,毕业于中国石油大学(北京)地质工程专业,现工作于中国石油新疆油田分公司勘探开发研究院,目前从事地震资料处理及储层预测等工作。。

摘要

针对现有常规去噪方法手段在参数选择及条件假设的不合理因素,导致去除的噪声剖面中残留部分有效信号(又称能量泄漏)的问题,本文提出了基于自适应K-SVD的能量泄漏恢复方法,实现残存有效信号的回收利用。首先采用小波阈值、f-x域去噪等传统方法处理,获得一个较高品质的样本进行字典的构造,再引入相关系数匹配准则对残差剖面进行稀疏分解,实现随机噪声与残留有效信号的二次分离。试验结果表明,本文算法能较好地兼顾地震信号的去噪与保幅,具有自适应性。

关键词： 能量泄漏 ; 优质样本 ; 相关系数匹配准则 ; 信噪分离

Abstract

In view of the problem of residual signal (also known as energy leakage) in the noise profile of the existing conventional denoising method in the parameter selection,this paper proposes the energy leakage recovery method of adaptive K-SVD to achieve effective recycling of signals.Firstly,a high-quality sample is obtained by using traditional methods such as wavelet threshold and fx domain denoising to perform dictionary construction,and then the correlation coefficient matching criterion is introduced to sparsely decompose the residual profile so as to realize the secondary separation of random noise and residual effective signal.The experimental results show that the proposed algorithm can better balance the denoising and amplitude preservation of seismic signals and has adaptability.

Keywords： energy leakage ; high-quality samples ; correlation coefficient matching criterion ; signal-to-noise separation

PDF (9142KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

何真, 曹思远, 郝婳婕, 段绪远, 叶力夏提, 张韵佩. 基于自适应K-SVD的能量泄漏恢复研究. 物探与化探[J], 2020, 44(2): 362-371 doi:10.11720/wtyht.2020.1229

HE Zhen, CAO Si-Yuan, HAO Hua-Jie, DUAN Xu-Yuan, , ZHANG Yun-Pei. Research on energy leakage recovery of adaptive K-SVD. Geophysical and Geochemical Exploration[J], 2020, 44(2): 362-371 doi:10.11720/wtyht.2020.1229

0 引言

在地震资料处理中,随机噪声的压制研究一直都是热点与难点,处理结果的好坏直接影响着后续解释与储层预测等工作。目前生产实际中常用的去噪手段有小波变换、奇异值分解(SVD)、F-X域去噪、Radon变换、三维块匹配算法(BM3D)等一系列方法。其中,小波变换^[1]作为一种经典的数学变换,对于不同特征的带噪信号采用软阈值函数或硬阈值函数进行降噪处理,虽有一定效果但仍然有待提高,该方法在提出以后也不断有专家学者对其进行讨论与改进^[2,3]。奇异值分解算法(SVD)^[4]借助地震信号之间的相关性与随机噪声之间的不相关来进行信噪分离,根据能量的不同分为若干奇异值,大奇异值对应有效信号,小奇异值对应噪声,但地震同相轴斜率较大时则不适用该方法。后来刘志鹏、赵伟等^[5]提出局部频率域SVD对地震波进行去噪,基于有效频带频谱存在较好的相关性来进行噪声压制。f-x域去噪^[6,7]方法假定相干信号在f-x域可预测而随机噪声不可预测的原理对信号和噪声加以区分来达到去噪的目的。Radon变换去噪主要有线性、抛物线和双曲线拉东变换3种去噪方法^[8,9,10],根据同相轴的不同时差特征来选取合适的变换类型进行去噪。三维块匹配算法^[11,12](BM3D)是近些年热度较高的图像去噪算法之一,它将噪声图像沿垂直与水平方向搜索相似块并分组,再进行协同滤波,该方法解决了图像去噪后模糊的问题,保留了图像的边缘细节,在一定程度上优于传统的空间域降噪技术和频率域降噪技术。

尽管去噪方法和理论在不断改进和完善,但在常规处理中由于参数选择的不合理或者基于一些不合理假设导致在对信号进行降噪的同时难免会“伤及”一部分有效信号,即残差剖面中存在有效信号的情形,称其为能量泄漏现象。对于噪声中含有的有效信号部分很少有人进行研究。本文提出基于优质样本的学习字典进行能量恢复研究,优质样本构建的字典包含较多的原地震信息以及较少的随机噪声,并且与残差剖面中残留的有效信号有较强的相关性,可以有效地将其提取并回收利用。

1 自适应K-SVD算法能量泄漏恢复原理

1.1 K-SVD算法

实际地震数据可表示为:

(1)

f = s + σ * n, (1)

其中:f为含噪地震数据;s为有效信号;n为随机噪声;σ为噪声强度,一般不可知。压制随机噪声的过程就是将σ*n这一部分从地震数据f中去除,现用任意手段将信号f进行去噪处理得到:

(2)

D_noise = Method (f),

其中:D_noise表示去噪后信号,Method表示任一种去噪方法,可以是小波阈值去噪、SVD分解、中值滤波等。

一次去噪后信号D_noise即本文研究的部分,将其写成如下形式:

(3)

D_noise = D * X,

式中:D_noise∈R(m×n),D∈R(m×k),X∈R(k×n),k是训练字典中系数的长度,D是超完备字典,X是对应D的稀疏稀疏矩阵,训练数据D_noise={d_noise_i} $_{i = 1}^{n}$ 表示为训练样本的集合,字典D={d_i} $_{i = 1}^{k}$ 的每一列表示一个原子,矩阵X={x_i} $_{i = 1}^{n}$ 为n个稀疏向量的集合。要找到一个D让X尽量稀疏,即:

(4)

\begin{array}{l} \min_{D . X} ‖ D_noise - DX ‖_{F}^{2} \forall i, ‖ x_{i} ‖_{0} \leq T_{0} \\ i = 1,2, \dots, n \end{array}

T₀为稀疏表示中非零分量个数的最大值。字典D的更新是逐列进行,假设字典D和稀疏矩阵X都是固定的,当更新字典D中第h列原子d_h时,则有:

(5)

\begin{array}{l} ‖ D_noise - DX ‖_{F}^{2} = ‖ D_noise - \overset{h}{\sum_{j = 1}} d_{j} x_{T}^{j} ‖_{F}^{2} \\ = ‖ (D_noise - \overset{h}{\sum_{j \neq h}} d_{j} x_{T}^{j}) - d_{h} x_{T}^{h} ‖_{F}^{2} \\ = ‖ E_{h} - d_{h} x_{T}^{h} ‖_{F}^{2} \end{array}

为了使误差矩阵E_h尽可能小,利用SVD分解计算可选的d_h与相应的系数 $x_{T}^{h}$ 能够有效最小化E_h的误差。定义ω_h存储用到 $x_{T}^{h}$ 的非零索引,即:

(6)

ω_{h} = {i | 1 \leq i \leq h, x_{T}^{h} (i) \neq 0},

再定义Ω_h矩阵,大小为n×|ω_h|,在(ω_h(i),i)取值为1,其余位置为0。令向量 $x_{R}^{h}$ = $X_{T}^{h}$ *Ω_h,D_nois $e_{R}^{h}$ =D_noise*Ω_h, $E_{R}^{h}$ =E_h*Ω_h。其中 $x_{R}^{h}$ 长度为|ω_h|,仅保留了系数中非零成分,D_nois $e_{R}^{h}$ 与 $E_{R}^{h}$ 大小均为n×|ω_h|,且只保留与 $X_{T}^{h}$ 中非零部分对应的列。此时上式可转化为:

(7)

‖ E_{h} Ω_{h} - d_{h} x_{T}^{h} Ω_{h} ‖_{F}^{2} = ‖ E_{R}^{h} - d_{h} x_{R}^{h} ‖_{F}^{2},

对 $E_{R}^{h}$ 进行SVD分解:

(8)

E_{R}^{h} = UΔ V^{T} 。

式中:U的列向量是E( $E_{R}^{h}$ )^T的特征向量,V的列向量是E^T $E_{R}^{H}$ 的特征向量。Δ是按照奇异值递减的顺序组成的对角矩阵:

(9)

Δ = [\begin{array}{l} Σ & 0 \\ 0 & 0 \end{array}],

(10)

Σ = diag (σ_{1}, σ_{2}, σ_{3}, \dots, σ_{r}) 。

σ_r(r=1,2,3,…,r)表示矩阵的所有非零奇异值,r为矩阵的秩。

将矩阵U第一列U₁替换字典D第h列,V₁表示矩阵V第一列,则系数矩阵 $x_{R}^{h}$ 可表示为V₁*Δ(1,1)(最大奇异值),此时字典d_h更新完成,将字典按照这样的方式不断迭代更新,直到结果满意为止。

经K-SVD处理构建得到一次去噪后剖面的字典信息,定义为H-D,再使用该字典对一次去噪后的残差剖面RES进行稀疏分解,分离得到残差剖面中残存的有效信号 $\bar{RES}$ ,即

(11)

\bar{RES} = KSVD (RES),

将提取残差剖面的有效信号部分 $\bar{RES}$ 补回去噪后结果D_noise即完成了处理:

(12)

Rec = D_noise + \bar{RES} 。

1.2 正交匹配追踪(OMP)

OMP算法是根据匹配度选择最优的原子进入原子集,它在更新信号残差时,减去的是当前残差在所有已被选择的原子上的正交投影,这样就使得该算法在迭代过程中避免对已选的原子的重复挑选。并且通过求解一个最小二乘问题找到原信号的最优近似解,再对残差进行更新迭代,最后用原子的线性组合来表征信号。具体算法如下:

1)输入m×n的字典D以及需要稀疏重构的信号y,给定预设误差ε,稀疏度K;

2)初始化残差ε₀=y,索引集Λ₀=Ø,迭代次数t=1;

3)将字典D每列与残差ε_t做内积,选择内积最大的一个即信号在原子垂直投影长度中最长的一个,并记录其索引ϕ_t;

4)令Λ_t=Λ_t_-1∪{ϕ_t};

5)求解最小二乘,得到 $\tilde{X}$ =argmin‖y-DX‖₂;

6)更新残差ε_t=y-D $\tilde{X}$ ,t=t+1;

7)判断是否满足ε_t<ε,若满足,停止迭代,否则返回3)。

1.3 相关系数匹配准则

一般基于字典学习的迭代停止准则采用残差信号能量的迭代停止标准^[13],但残差剖面中随机噪声成分较多,有效信号相对较少,该准则过早终止了OMP迭代,导致有效信号提取不充分。本文引用一种新的基于相关系数匹配的停止准则^[14]。该准则定义匹配因子 $M_{γ}^{t}$ 是第t次迭代以后二次残差信号(残差剖面的残差)矩阵RRES^t^-1和字典原子d_γ之间的相关系数,即:

(13)

M_{γ}^{t} = \frac{< RRE S^{t - 1}, d_{γ} >}{‖ RRE S^{t - 1} ‖_{2}} 。

匹配因子的最大值会随着信号的不断迭代而减小,设定一个阈值T_m,假定匹配因子 $M_{γ}^{t}$ 最大值小于阈值T_m时,残差剖面中的有效信号基本被提取完毕,即:

(14)

\max_{d_{γ} \in D_{t}} | M_{γ}^{t} |^{2} \leq T_{m} 。

此时停止OMP迭代,输出结果。

1.4 本文算法实现流程

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 本文算法流程图

Fig.1 Algorithm flow chart

2 实验结果与分析

实验分为对模型资料的处理以及实际资料的处理,首先通过楔形模型和倾斜模型的处理证明本文方法的正确性和有效性,再对实际资料进行分析、去噪、能量恢复获得高信噪比高保真度的地震剖面。

实验中为统一标准,设置地震数据块大小为8×8,字典中的原子数为256,初始化字典采用DCT字典。在评价去噪效果时采用能量叠加法^[15]计算信噪比,公式为:

(15)

E_{S} = N \overset{M}{\sum_{i = 1}} s_{i}^{2} = N \overset{M}{\sum_{i = 1}} {[\frac{\overset{N}{\sum_{j = 1}} x_{ij}}{N}]}^{2} = \frac{1}{N} \overset{M}{\sum_{i = 1}} {[\overset{N}{\sum_{j = 1}} x_{ij}]}^{2},

(16)

SNR = \frac{E_{S}}{E_{N}} = \frac{E_{S}}{E - E_{S}} = \frac{\overset{M}{\sum_{i = 1}} {[\overset{N}{\sum_{j = 1}} x_{ij}]}^{2}}{N \overset{M}{\sum_{i = 1}} \overset{N}{\sum_{j = 1}} x_{ij}^{2} - \overset{M}{\sum_{i = 1}} {[\overset{N}{\sum_{j = 1}} x_{ij}]}^{2}} 。

式中:M、N分别为采样点数和地震道数,SNR为信噪比,E为地震道总能量,E_N为噪声总能量,E_S为有效信号能量。

2.1 楔形模型去噪分析

图2a是楔形模型,道数30道,采样点数500,子波主频35 Hz;图2b是加入了50%的随机噪声,其信噪比为2.18 dB。现在利用KSVD算法分别对无噪和带噪情况下的模型求取特征字典,得到图2c(带噪字典),图2d(无噪字典)。不难发现在无噪情况下,特征字典提取的都是有效信号的成分,而噪声情况下的字典则较为“模糊”。现对加噪信号分别采用带噪和无噪的字典进行稀疏求解,得到去噪结果图2e和图2f。信噪比分别为18.34 dB和20.05 dB,两种方法对噪声均有较好的压制作用,但对比残差剖面图2g、h和单道15道信号图2i、j、k、l来看,采用无噪字典去噪的剖面在信号的去噪效果和保幅性上更胜一筹。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 楔形模型试验结果

a—楔形模型;b—加噪楔形模型;c—噪声情况下构建的字典;d—无噪声情况下构建的字典;e—噪声字典去噪剖面;f—无噪声字典去噪剖面;g—噪声字典去噪后残差;h—无噪声字典去噪后残差;i—楔形模型第15道信号;j—加噪楔形模型第15道;k—噪声字典去噪后第15道;l—无噪字典去噪后第15道

Fig.2 Wedge model test results

a—wedge model;b—noise wedge model;c—dictionary constructed under noise;d—noise dictionary denoising profile;e—noise dictionary denoising profile;f—noiseless dictionary denoising profile;g—noise dictionary denoising residuals;h—noiseless dictionary denoising residuals;i—wedge model 15^th signal;j—noise wedge model 15^th;k—noise dictionary denoising after 15^th;l—No15 after noiseless dictionary denoising

2.2 倾斜模型能量泄漏恢复效果分析

上文2.1中通过楔形模型证明利用一个优质的学习样本(即无噪声或噪声尽可能少的情况下的剖面)进行字典学习对噪声压制以及信号重构的效果要优于传统对带噪信号进行学习(即普通样本)的效果。接下来采用倾斜模型,对其加入一定比例随机噪声,先利用传统去噪手段对其进行去噪处理,再分别对去噪后剖面和原始带噪剖面作为样本进行字典学习,将学习到的字典对残差剖面进行稀疏分解,实现信号和噪声的二次分离,最后将分离出的有效信号补回去噪后的剖面并分析讨论两种方法的恢复结果。

图3a为倾斜模型,道数81道,采样点数801,子波主频55 Hz;图3b加入65%的随机噪声,信噪比为3.31 dB。利用传统小波变换对其进行去噪处理,得到图3c,信噪比为19.62dB,去噪后剖面信噪比有了大幅度的提升,但从残差剖面观察(图3d),去除的噪声中存在大量有效信号。以普通样本和优质样本为基础进行字典学习对残差剖面进行处理,得到图3e和图3f的信噪分离结果。从分离结果可以看出,以优质样本学习的字典在残差剖面中信噪分离的效果更加彻底,将其分别补回去噪后剖面得到能量恢复结果(图3g、h),信噪比分别为18.13dB和21.25 dB。可见本文提出的方法对有效信号的提取更加精准和可靠,下文的单道信号分析中更能清晰地体现这一结论。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 倾斜模型剖面试验结果

a—倾斜模型;b—加噪倾斜模型;c—小波阈值去噪剖面;d—小波阈值去噪残差;e—普通样本分离的有效信号;f—优质样本分离的有效信号;g—普通样本能量恢复结果;h—优势样本能量恢复结果

Fig.3 Tilt model section test results

a—tilt model;b—charging tilt model;c—wavelet threshold denoising profile;d—wavelet threshold denoising residual;e—effective signal for normal sample separation;f—effective signal for good sample separation;g—ordinary sample energy recovery results;h—high quality sample energy recovery results

对图3倾斜模型时间域第75道,观察普通样本的字典学习恢复效果(图4d)和优质样本的字典学习恢复效果(图4e),优质样本学习下恢复的结果形态更加接近原始信号,而普通样本学习下恢复的结果中含有较多的随机噪声,信噪分离程度相对不够理想。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 倾斜模型第75道信号试验结果

a—倾斜模型;b—倾斜模型加噪;c—小波阈值去噪结果;d—普通样本能量恢复结果;e—优质样本能量恢复结果

Fig.4 Test results of the 75^th signal of tilt model

a—tilt model;b—tilt model noise;c—wavelet threshold denoising result;d—normal sample energy recovery result;e—good sample energy recovery result

为更好地体现本文方法的实用性,现将带噪倾斜剖面分别采用f-x域去噪、改进小波阈值^[16]、ASVD^[17]3种方法对其进行一次去噪,再分别对3个去噪剖面进行字典构造进而提取残差剖面的有效信号,并对结果进行分析。

由于补偿恢复后的结果从时间域剖面中很难看出差异,所以对以上3种恢复结果只从单道信号进行观察讨论。从去噪结果看,3种方法均有较好的效果,噪声基本被去除,但f-x域去噪相对于改进小波阈值和ASVD方法,在去除噪声的同时也伤及了大量有效信号。仍然以恢复后的第75道为例,见图5g、h、i,对f-x域去噪、改进小波阈值以及ASVD三种去噪结果进行能量恢复后随机噪声含量均比较小,但在有效信号振幅能量上后两者恢复得更好,更加接近原始模型资料。经对比也不难发现:原始剖面降噪效果越好,有效信号保留越多,在构建字典原子时信号的特征越能被更好地提取和表达,从而实现更加理想的恢复效果。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 本文方法实用性试验结果

a—f-x域去噪剖面;b—f-x域去噪残差;c—改进小波阈值去噪剖面;d—改进小波阈值去噪残差;e—ASVD去噪剖面;f—ASVD去噪残差;g—f-x域去噪后能量恢复结果(第75道);h—改进小波阈值去噪后能量恢复结果(第75道);i—ASVD去噪后能量恢复结果(第75道)

Fig.5 Result of the practical test of the method

a—f-x domain denoising profile;b—f-x domain denoising residual;c—improved wavelet threshold denoising profile;d—improved wavelet threshold denoising residual;e—ASVD denoising proflie;f—ASVD denoising residual;g—f-x domain denoising energy recovery result(75^th channel);h—improved wavelet threshold denoising energy recovery result(75^th channel);i—ASVD denoising energy recovery result(75^th channel)

2.3 实际资料处理及效果分析

选取某实际二维叠加数据进行处理分析。图6a为某区块实际地震资料,剖面上随机噪声较为发育,能量较弱的信号基本被湮没,给接下来的小层解释、储层预测等工作带来了较大难度。实际资料的一次去噪选取fx域和ASVD法,进行f-x域去噪,得到去噪结果图6b和残差剖面图6c,采用本文方法对残差剖面处理实现噪声信号二次分离,将分离的有效信号加回去噪后剖面得到能量恢复结果(图6d),同样地,采用上述方法得到ASVD的处理结果(图6e、f、g),对于一次去噪结果(图6b和图6e),两种方法均去除大部分随机噪声,使得有效波变得清晰,但残差剖面中仍能发现较多的有效信号。而恢复后的剖面图6d、g中同相轴信息被较好地提取(见红圈及红框部分),剖面变得更加自然、连续。不难发现,ASVD处理后再进行能量恢复的剖面在细节表现上要优于f-x去噪后再恢复的结果。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6 实际资料处理

a—某地区二维实际资料;b—f-x域去噪剖面;c—f-x域去噪残差;d—f-x域去噪后恢复剖面;e—ASVD去噪剖面;f—ASVD去噪残差;g—ASVD去噪后恢复剖面

Fig.6 Actual data processing

a—2D actual data in a region;b—f-x domain denoising profile;c—f-x domain denoising residual;d—f-x domain denoising recovery profile;e—ASVD denoising proflie;f—ASVD denoising residual;g—ASVD domain denoising recovery profile

3 结论与建议

本文基于K-SVD字典学习方法进行地震能量泄露恢复研究,通过楔形、倾斜模型以及实际资料处理分析得到以下结论:

1)以无噪或者噪声尽可能小的剖面作为样本进行字典学习、稀疏表示得到的去噪结果优于传统的字典学习方法。

2)利用一次去噪后的剖面作为样本进行学习,原始剖面的降噪性、保真性越好,在构建字典原子时信号的特征越能被更好的提取和表达,从而实现更加理想的恢复效果。

3)根据不同的地震信号可以构建不同的特征字典,具有较好的自适应性,并且在适用范围上有一定应用前景。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Donoho D

Denoising by soft-thresholding

[J]. IEEE Trans on Information Theory, 1995,41(3):613-627.