基于局部峰值约束的快速正交匹配追踪地震数据分解方法

杜泽源, 杨 森, 陶永慧, 蔡杰雄, 何兵红

(1. 中国石油化工股份有限公司 石油物探技术研究院，南京 211103；2. 中国石油大学(华东) 地球科学与技术学院，青岛 266580)

0 引言

地震数据包含了丰富的地下介质信息，通过对地震数据分解可以得到多个含有不同信息的地震资料，有利于提高对地下不同介质情况及数据特征的认识与辨别。地震数据分解是分析地震数据结构、属性特征，提取地质信息的有效途径。在地震信号处理领域，地震数据稀疏分解可以用于数据去噪、数据重建[1]以及时频分析[2]和地震反演等问题[3]。数据分解可以通过多种途径实现，傅里叶变换将时间域信号变换到频率域，使人们可以更好地分析信号的特点，在频率域对信号进行处理在很多情况下可以带来很多便利[4]。但是傅里叶变换只能对信号整体进行分析，不能获得局部的信息[5]，这对于非平稳地震信号的处理是不利的。小波变换以其多分辨率的特点能够得到更好的信号分解，小波分析的思想被广泛应用到地震数据分析中[6-7]。Mallat[8-9]提出了基于塔式多分辨率分析的小波分解与信号重构的快速算法(Mallat算法)；许康生[10]研究了基于小波变换的地震信号多尺度分解，证明了重建地震信号无畸变和虚假信号产生，而且多尺度分解可用于地震记录的弱震相时间信息的提取；陈飞等[11]利用小波分析进行地震信号多尺度分解及去噪处理，通过数据测试证明能够得到较好的去噪效果；杜丽英[12]提出小波包最佳基地震数据分解与重建，过最佳小波基的选择，使在空间域局部性很好的情况下，频域的局部性也达到最佳状态；周静毅等[13]研究了基于小波变换的地震奇异性分析，利用小波变换良好的时间-频率局部化特征，通过地震数据奇异性属性分析来划分地层旋回。但小波变换母小波的选择以及分解的层数，对信号多尺度分解的结果影响较大，有时候难以获得理想的尺度分解数据。

稀疏表示理论为信号分析处理带来了一次革命性的推动，它的核心是通过稀疏分解，借助字典中少量原子的线性组合形式对原信号所涵盖的信息进行表达。Mallat等[14]提出使用过完备原子库用于信号稀疏分解的匹配追踪算法。匹配追踪算法具有自适应选择匹配原子的优点[15-16]，在地震信号处理领域的应用发展迅速。何胡军等[17]研究了基于匹配追踪算法子波分解技术识别薄互层储层，该技术能较好地识别薄互层储层，分辨出储层在空间的展布规律；张繁昌等[18]利用基于匹配追踪的地震瞬时谱分解方法实现了三角洲砂岩尖灭先的识别。匹配追踪是一种贪婪算法，在信号分解过程中容易出现重复选择原子的情况，影响计算效率，许多学者也提出了改进算法。Wang[19]提出了多道匹配追踪算法，并应用于地震数据的分解；汪瑞良等[20]研究了快速动态匹配追踪算法，实现了匹配追踪高分辨率时频谱计算。正交匹配追踪(Orthogonal Matching Pursuit，OMP)算法在每次迭代计算过程中保证选择的原子相互正交，计算效率更高[21]。快速有效地进行匹配追踪计算是一个挑战性问题，对于最初的正交匹配追踪算法，由于原子字典规模很大，对原子遍历搜索，计算成本依然很高，对于正交匹配追踪算法的优化依然是研究的热点。

针对正交匹配追踪在通过内积实现子波与信号匹配时不能保证子波与信号具有相同的位置及振幅的问题，避免在子波匹配时对字典内所有子波都进行扫描而大大降低计算效率。笔者采用局部峰值约束的正交匹配追踪方法来解决这一问题，构建信号与子波峰值距离的目标泛函，通过快速扫描算法进行高效的求解计算。模型测试结果表明，该方法可以高效实现地震数据的高精度分解与重构，验证了本文算法的有效性和高效性。

1 理论方法

1.1 正交匹配追踪求解方法原理

匹配追踪是基于稀疏表示理论发展起来的，稀疏表示为地震信号处理及地震反演等方面带来了很多便利。地震信号的稀疏表示可以通过几个原子(基本波形)的线性加权组合表示：

(1)

其中：s(t)为地震信号；ω(t)为信号原子(基本子波)；ωk(t)是从一个原子库(子波字典)中选取的；原子库Φ=[ω1(t),ω2(t),…,ωK(t)]是由信号原子为列向量构成的维度为N×K的矩阵；N为子波长度；αk为信号原子对应的系数；α=[α1,α2,…,αK]为信号在原子库Φ中的表示系数。

匹配追踪算法的主要思想是，首先建立一个过完备原子库，在每次迭代求解过程中，从原子库中选出与信号最为匹配的原子，来稀疏逼近信号，进而求出信号与所选原子的残差；然后从原子库中选出与残差信号最匹配的原子，同时计算新的残差；重复多次迭代计算以上过程，直到残差小于设定的阈值或迭代次数达到设定的上限迭代次数时，将选出的最佳原子进行线性组合就可以表示原始信号。但是匹配追踪存在一个问题，在匹配分解过程中很可能会重复选择原子，这严重影响计算效率。正交匹配追踪算法要求每次求得的信号残差与之前选出的所有原子相正交，这样在匹配时就会避免重复地选择原子，能够更快地收敛到全局极值[21-22]。

所述正交匹配追踪算法的目标泛函可表示为：

(2)

其中，σ为迭代终止阈值。

在稀疏求解问题中，如果原子ωk(t)正交，系数αk可以用信号s(t)与原子ωk(t)的内积来计算得到，即：

(3)

匹配迭代的过程为：首先从匹配原子库Φ中选出与给定地震信号s(t)匹配最佳的原子ω1(t)，匹配最佳的要求则是它们的内积[s(t),ω1(t)]为匹配原子库中所有原子与该地震信号s(t)内积最大的一个：

|[s(t),ω1(t)]|>|[s(t),ωk(t)]|,k≠1

(4)

即，第一次迭代后，信号s(t)被分解为最佳原子ω1(t)上的分量和信号残差R1(t)，

[s(t),ω1(t)]ω1(t)+R1(t)

(5)

然后对残差信号R1(t)重复上述相同的匹配分解过程：从Φ中筛选出与R1(t)最为匹配的第二个原子，即第二次迭代可表示为：

[R1(t),ω2(t)]ω2(t)+R2(t)

(6)

不断重复上述过程，第n次分解后，获得残差信号Rn(t)：

[Rn-1(t),ωn(t)]ωn(t)+Rn(t)

(7)

所以经过n次分解后，信号可表示为式(8)。

(8)

其中R0(t)为原始地震信号s(t)。式(8)表达的意思即为经过n次迭代分解计算后，原始信号s(t)可以用n个匹配原子与残差Rn(t)的合成来表示。

由于残差‖Rn(t)‖的衰减特性，对有限长度的信号，‖Rn(t)‖随着迭代次数n的增大而衰减至小于阈值σ，即信号的稀疏分解跟信号的近似程度足够好或残差信号满足设定的阈值时，迭代停止。一般，用少数M个原子就可以表示信号主要成分，满足数据处理的需求，且阈值σ是很小的数，则地震信号(可稀疏表示)为最佳原子的线性组合式(9)。

(9)

1.2 子波字典构建

匹配追踪算法能够快速有效计算的关键之一是，选取合适的信号原子构建原子库。经典的匹配追踪算法是基于Garbor原子来构建过完备匹配原子库。为适应地震信号不同的结构特征，对于原子库可以由一系列能够表达地震信号特征的子波波形构成，即过完备匹配子波库是一系列匹配子波构成的子波矩阵。

对于地震数据处理，子波字典实质上就是一系列基本子波ωδ(t)进行适当的时移、调制和相位变化得到的子波矩阵。若ωδ(t)为子波库中满足匹配条件的一个子波，其控制参数用δ(τ,f,φ,c)来表示，子波函数记为：

ωδ(t)=ω(t-τ)exp{i[πf(t-τ)/c+φ]}

(10)

其中：τ表示时移；f表示中心频率；φ表示相位；c表示尺度因子；i表示虚数单位。

子波字典的构建有很多子波可选择，可以根据所分析信号的特征及处理目标进行选择。对于地震数据来说，可以选择Ricker子波、高斯函数或Morlet小波等基本子波构建子波字典。Ricker子波是比较常用的子波，由于其波形简单，延迟时间短，时频分辨率较高，控制参数改变灵活，易于构建波形丰富的子波字典，充分表达地震数据[23]。另外，Morlet小波有较好的尺度变化能力，在时间和频率的局部化之间有着很好的平衡，切与地震信号有较好的相似性，也可以选取Morlet小波创建子波字典[24]。

1.3 局部峰值约束正交匹配追踪

正交匹配追踪的目标是最大化所选匹配子波与信号残差之间的内积，然而在内积实现时往往只是在全局时间尺度上度量原子与残差之间的信号相似性，并不能保证原子波峰与残差波峰具有相同的位置。

为了解决这一问题，采用局部峰值约束的方法来确定匹配子波参数，将最小化信号残差波峰与原子峰值距离为目标泛函。为了更有效地对信号及原子波峰进行扫描，首先扫描信号包络，获取其局部极大值来限定信号的扫描范围，然后再寻找信号峰值并扫描原子进行匹配。因此定义了以下目标函数：

(11)

具体地，对信号进行Hilbert变换，构造复数信号，通过时间序列确定信号采样率和信号长度，利用复数信号求取信号的瞬时包络与瞬时相位，则τi可粗略估计为瞬时包络局部极大值峰值处所对应的时间ti。在复数域内对子波库中的原子进行频率参数的扫描，对于某次迭代确定出的一组中心时间τb，得到信号在这些时间处的频率先验值fb，假设基本子波为零相位子波，则可利用τb和fb确定若干的先验信息点(τb,fb,0)。以这些先验信息点为中心，选取频率扫描范围，同时搜索多个原子，来提高匹配追踪效率，确定匹配子波的波峰频率fi。

在确定了最佳匹配子波ω(t)的中心时间和主频后，需要确定这些匹配子波的振幅和相位。根据残差能量尽可能小的原则，可以利用残差能量阻尼最小二乘方法[25]计算这些匹配子波的振幅和相位信息。

1.4 地震数据多尺度分解

地震信号稀疏处理一般可以表示为线性形式[1]：

dobs(t)=Fs(t)+o(t)

(12)

其中:dobs(t)为观测地震数据；o(t)为噪声；F为有界线性算子，对于不同的数据处理问题有不同的表达形式。对于本研究中，F表示多尺度分解算子；s(t)为稀疏求解数据。

地震资料可以分解成不同的分量(或尺度)。重建的地震信号可以表示为地震信号分量的线性叠加：

s(t)=sδ(f1)(t)+sδ(f2)(t)+…+sδ(fC)(t)

(13)

其中，δ(fi)代表不同的频带尺度。子波字典中子波由δ(τ,f,φ,c)参数控制包含频率因子，将频带尺度的子波归为一类，则子波字典可以分为相应频率尺度下Φδ(fi)。

地震信号稀疏表示的解为s(t)=Φα，则稀疏约束反演问题在相应频带尺度下可表示为求解以下泛函：

(14)

对地震数据稀疏分解则可以直接提取不同的频率分量，利用正交匹配追踪算法求解多尺度分解目标泛函(14)，即可将地震数据分解为不同频率尺度下的分量，基于本方法的数据分解有较好的时频分辨率，能够较好地服务于多尺度数据分析或多尺度反演。

2 模型测试

为验证文中正交匹配追踪算法的效果，设计几个基本波形叠加得到单道信号，如图1(a)所示。然后利用正交匹配追踪信号分解算法对合成地震信号进行分解，将分解后的信号进行重构，得到重构地震信号，见图1(b)。

对合成信号速度分解求解迭代次数为23次，从图1可以看出,重构信号与原始合成信号波形基本一致，图1(c)为原始信号与重构信号的残差，信号残差小于2‰，满足计算精度要求。图2为将原始信号与重构信号匹配对比，也可以看出原始信号与重构信号基本吻合，由此说明本文分解方法能够有效地匹配出最佳子波，并精确计算子波的振幅及相位。通过模型测试验证了本文稀疏分解重构算法能够将信号精确分解重构，效果较好。

图1 合成信号分解重构的结果Fig.1 The result of decomposition and reconstruction of synthetic signal(a)合成信号;(b)重构信号;(c)信号残差

图2 原始信号与重构信号对比Fig.2 Comparison of original signal and reconstructed signal

匹配追踪有较好的信号分解与重构精度，同时具有较高的时频分辨率。基于匹配追踪的短时傅里叶变换可以对每一个匹配子波进行变换，时窗调节比较灵活，相对于短时傅里叶变换时窗固定的特点，匹配追踪算法的时频分辨率更高。Wigner-Ville分布虽然相对于短时傅里叶变换时频分辨率提高，但是易产生交叉干扰项。分别用短时傅里叶变换、Wigner-Ville分布及基于匹配追踪的方法分析信号的时频谱，结果如图3所示。

图3 不同算法下的信号时频谱Fig.3 Time spectrum of signals under different algorithms(a)短时傅里叶时频谱;(b)Wigner-Ville时频谱; (c)基于OMP的短时傅里叶时频谱;(d)基于OMP的Wigner-Ville时频谱

从图3可以看出,基于OMP的短时傅里叶变换时频谱(图3(c))相对于常规短时傅里叶变换时频谱(图3(a))分辨率更高。常规Wigner-Ville分布时频谱(图3(b))比常规短时傅里叶变换时频谱分辨率高，但是产生了交叉项，而通过正交匹配追踪算法改进的Wigner-Ville分布时频谱(图3(d))消除了干扰项，同时保留了较高的时频谱分辨率，能量更加聚焦。通过对比说明基于匹配追踪的时频谱具有较高的时频分辨率和能量聚焦特性，能够准确表示地震信号的局部频谱特征。

3 实际资料分析

采用某实际资料分析本文方法应用效果，选取一条二维剖面(图4)。地震数据共131道，每道数据为0.5 s，采样间隔为2 ms。使用本文正交匹配追踪数据分解方法对图4中地震数据进行分解与重构，图5展示了原始地震剖面(图5(a))与重构后地震剖面(图5(b))，通过对比图5(a)和图5(b)可以看出，重构后的剖面与原始剖面基本完全一致，基本恢复了原始数据中包含的细节部分。

图4 原始地震剖面Fig.4 Original seismic section

图5 原始地震剖面与OMP分解重构剖面Fig.5 Original seismic section and OMP decomposition reconstruction section(a)原始地震剖面；(b)OMP分解重构剖面

抽取地震剖面中第75道数据，如图6(a)所示，在0.33 s处有一含气储层。对地震道进行时频分析，图6(b)为其时频谱，可以看出在含气储层位置处谱能量较为聚焦，时频分辨率较高，可以清楚地区分上下地层。通过时频谱可以非常直观的了解含气储层地震响应的频带特征。

图6 第75道地震数据及其时频谱Fig.6 The 75th trace seismic data and its time-frequency spectrum(a)第75道地震数据；(b)地震道时频谱

使用本文多尺度分解方法，对地震剖面进行尺度分解。由图7可以看出，分解后的地震剖面符合原始剖面中的地质特征，同相轴与原始剖面中形态一致，没有假频产生，频谱的混叠也能有效分离，也可以看出不同的地层有不同的频带相应特征。在第75道位置处有一口井，声波测井数据见图7，分解后的地震同相轴指示地层与测井速度分层相吻合。因此，可以不同地层地震响应频带特征为依据，使用分解后的地震数据在目标优势尺度下进行数据分析或地震反演，会得到更好的效果。通过地震数据分解，可以为面向各种目标的数据处理或地震反演问题选择提供可靠的数据分量。

图7 地震数据多尺度分解剖面Fig.7 Multi-scale decomposition sections of seismic data(a)主频为10 Hz的单频剖面；(b)主频为16 Hz的单频剖面； (c)主频为22 Hz的单频剖面；(d)主频为28 Hz的单频剖面

4 结论

匹配追踪算法是一种贪婪算法，难以避免原子的重复选择，正交匹配追踪算法要求每次迭代后的残差信号和之前选出的所有原子相正交，这样在匹配时就会避免重复地选择原子。但是正交匹配追踪在内积的实现时，往往只是在全局时间尺度上度量原子与残差之间的信号相似性，并不能保证原子波峰与残差波峰具有相同的位置，对信号原子的扫描也存在冗余。针对这一问题，笔者研究了基于局部峰值约束的正交匹配追踪方法，并在原子扫描时进行优化，避免过多的选取原子，影响效率。在算法实现的基础上，通过模型数据与实际资料测试证明了本文数据分解与重构方法的有效性。

地震数据稀疏表示分解方法能够高效精确地将地震数据分解为不同频带尺度，在数据稀疏分解时够克服频谱混叠和波的干涉效应，提高时频分辨率，并且不会引入假频。得到的多尺度(分量)地震数据，可以为数据分析或多尺度地震反演提供可靠的数据支持。