0 引言

随着油气勘探程度的不断提高，勘探目标越来越复杂，对地震资料分辨率的要求也越来越高。近年来，随着大数据和人工智能技术的发展，机器学习也与地球物理领域有机结合，提高了传统地震信号处理的运算效率，但在分辨率方面的研究和应用还很少。

匹配追踪(Matching Pursuit，MP)算法^[1-3]是一种用于信号处理和机器学习的技术，特别是在稀疏表示和信号重构方面，能够把地震信号分解为一系列小波集合，实现信号的时频分解。但是，当前使用匹配追踪算法分解和重建多子波地震信号尚存在两方面问题^[4]：①地震子波旁瓣振荡严重导致时域追踪精度不足，地震数据分辨率降低；②传统匹配追踪算法分解地震道后，频谱成分众多、数据冗余，不能针对性地凸显目标地层的特征，缺乏智能化重建方法。

1901年，Pearson^[5]首次提出了主成分分析(Principal Component Analysis，PCA)方法，描述了一种通过维度变换把多元问题转化为更少元素的数学统计方法，用于提取数据特征和建立数理模型。但是PCA是基于线性变换的降维方法，仅适用于线性数据。为了更好地处理非线性数据，Schölkopf等^[6]将PCA方法扩展到非线性领域，并引入核函数进行映射，形成了核主成分分析(Kernel Principal Component Analysis，KPCA)方法。KPCA方法更适合处理非线性数据^[7]，通过非线性映射把原始空间的数据投影到高维空间，然后在高维空间提取数据特征，再用PCA方法进行数据处理。

本文将匹配追踪算法与KPCA方法应用于地震资料处理，提出了地震信号高分辨率处理方法。首先利用匹配追踪算法通过稀疏分解不断迭代，得到地震信号的最有效信息；然后，将子波替换为宽带Ricker子波进行整形处理，有效压制子波旁瓣，提高地震剖面分辨率；最后，采用KPCA方法将原始地震信号非线性映射到高维空间，在高维空间重建地震信号，消除冗余信息，可更准确地识别薄互层。最终实现地震信号的高精度、高效率、智能化重建，在保持原有信号准确性的同时，大大提高资料分辨率，为断层识别、层位追踪、薄层刻画等后续应用提供可靠的技术支持。

1 原理及实现方法 1.1 匹配追踪算法

在多子波地震道模型中，原始地震信号能够转化为各种相关联的子波原子集合，可以用以下形式刻画

$ S\left(t\right)={R}_{i}\left(t\right)\sum\limits_{i=1}^{I}{W}_{i}\left(t\right)+N\left(t\right) $

(1)

式中：$ S\left(t\right) $为地震信号；t为时间；$ {W}_{i}\left(t\right) $表示第$ i $个反射层的地震子波；$ {R}_{i}\left(t\right) $表示第$ i $个反射层的反射系数；$ N\left(t\right) $表示噪声项；I表示反射层总数。

匹配追踪算法从大量的Morlet小波字典中，自适应地筛选一系列与原始信号匹配的原子^[8]。Morlet小波能够精确地刻画实际地震道中存在的高频能量衰减现象^[9]，其时间域表达式为

$ \begin{align} m\left(t\right)=&\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{\mathrm{l}\mathrm{n}2}{{\mathrm{\pi }}^{2}}\frac{{\omega }_{{}^{\mathrm{m}}}^{2}{(t-u)}^{2}}{{\sigma }^{2}}\right)\times \\& \mathrm{e}\mathrm{x}\mathrm{p}\left[\mathrm{i}\left({\omega }_{\mathrm{m}}\right(t-u)+\psi )\right]\end{align} $

(2)

式中：$ {\omega }_{\mathrm{m}} $为平均角频率；$ u $是时移；$ \sigma $是尺度参数，决定波形形状；$ \psi $是相移。

小波的频谱可以表示为

$ \begin{array}{l}M\left(\omega \right)=\frac{1}{2}\sqrt{\frac{\mathrm{\pi }}{\mathrm{l}\mathrm{n}2}}\times \mathrm{e}\mathrm{x}\mathrm{p}\left[-\mathrm{i}(\omega u-\psi )\right]\times \\ \frac{\sigma }{{\omega }_{\mathrm{m}}}\mathrm{e}\mathrm{x}\mathrm{p}\left[-\frac{{\mathrm{\pi }}^{2}}{4\mathrm{l}\mathrm{n}2}\frac{{\sigma }^{2}(\omega -{\omega }_{\mathrm{m}}{)}^{2}}{{\omega }_{\mathrm{m}}^{2}}\right]\end{array} $

(3)

基于Morlet小波的匹配追踪算法步骤包括：①对输入地震信号进行Hilbert变换，构建复地震道；②刻画Morlet小波，参数$ u $为复地震道最大振幅包络处对应的时间，中心频率$ \omega $以及相移$ \psi $都可以用对应时刻的瞬时频率和瞬时相移代替^[10]；③在小波原子库里追踪最佳小波，计算小波振幅值；④从当前的残差信号中减去上一步输出的最佳小波，得到新的残差；⑤重复上述步骤，直至满足设定的迭代终止条件。

匹配追踪算法通过稀疏分解，利用小波原子库最大程度地逼近信号的有效分量，从而达到忽略噪声的效果。经过匹配追踪后，地震数据虽然保持了原始地震剖面的反射特征，但并不会改善原本地震子波旁瓣存在的振荡问题，导致时频追踪精度不足^[11]、数据的分辨率降低。因此，本文利用子波替换的方法对分解后的子波进行整形。

1.2 子波整形

地震资料的垂向分辨率可以通过子波主旁瓣比评价。当旁瓣能量被压制，地震资料的变化体现在同相轴减少，分辨率提高^[12]。

地震子波和反射系数序列在复赛谱上具有可分离性^[13]，首先将地震记录在频率域表示出来，可以用与式（1）类似的褶积模型描述^[14]

$ S\left(t\right)=W\left(t\right)\times R\left(t\right) $

(4)

式中：$ W\left(t\right) $为地震子波；$ R\left(t\right) $为地震反射系数序列。

根据对分辨率和保真度的需要构建均衡Ricker子波，再对地震信号进行反褶积处理，提取反射系数序列，把子波替换为提前构建好的宽带Ricker子波。宽带Ricker子波表示为

$ y\left(t\right)=\frac{1}{q-h}\left\{q\ \mathrm{e}\mathrm{x}\mathrm{p}[-{\left(\mathrm{\pi }qt\right)}^{2}]-h\ \mathrm{e}\mathrm{x}\mathrm{p}[-{\left(\mathrm{\pi }ht\right)}^{2}]\right\} $

(5)

式中$ q $和$ h $为Ricker子波频率参数的积分限。

相较于Morlet子波，宽带Ricker子波在提高地震信号的垂向分辨率方面具有显著的优势，在相同条件下，具有主瓣窄、旁瓣弱的特性^[15]，有效提高了地震信号的垂向分辨率，解决了子波旁瓣振荡严重导致的精度不足问题。两种子波的波形如图 1。

1.3 KPCA方法

KPCA和PCA有着相似的理论基础，都可以通过维度变换把多元问题转化为更少元素。KPCA是PCA在非线性领域的延伸。

1.3.1 PCA方法原理

PCA的基本思想是提取原始线性数据主要的特征代替原有数据^[16-17]，这些特征必须包含原始数据中绝大部分信息，并且彼此之间互不相关。假设有包含了$ n $个样本的样本集$ \{{\boldsymbol{x}}_{1}, {\boldsymbol{x}}_{2}, \cdots , {\boldsymbol{x}}_{n}\} $，其中每个样本点本身还具备$ p $个变量，则

$ \boldsymbol{X}=\left[\begin{array}{c}{{x}_{11}}&{{x}_{12}}&{\cdots }&{{x}_{1p}}\\ {{x}_{21}}&{{x}_{22}}&{\cdots }&{{x}_{2p}}\\ {}&{}&{ \vdots }&{}\\ {{x}_{n1}}&{{x}_{n2}}&{\cdots }&{{x}_{np}}\end{array}\right]=\left[{\boldsymbol{x}}_{1}\ {\boldsymbol{x}}_{2}\ \cdots \ {\boldsymbol{x}}_{p}\right] $

(6)

式中

$ {\boldsymbol{x}}_{i}=({x}_{1i}, {x}_{2i}, \cdots , {x}_{ni}{)}^{\mathrm{T}}, i=\mathrm{1, 2}, \cdots , p $

(7)

PCA的过程是将原始数据的$ p $个变量重新进行线性组合，从而得到一组$ p $个变量的新组合^[18]，即为

$ \boldsymbol{Z}=\left[\begin{array}{c}{{a}_{11}}&{{a}_{21}}&{\cdots}&{ {a}_{p1}}\\{ {a}_{12}}&{{a}_{22}}&{\cdots}&{ {a}_{p2}}\\{ }&{}&{ \vdots}&{ }\\{ {a}_{1p}}&{{a}_{2p}}&{\cdots }&{{a}_{pp}}\end{array}\right]\left[{\boldsymbol{x}}_{1}\ {\boldsymbol{x}}_{2}\ \cdots \ {\boldsymbol{x}}_{p}\right]={\boldsymbol{Q}}^{\mathrm{T}}\boldsymbol{X} $

(8)

其中

$ \boldsymbol{Q}=\left[\begin{array}{c}{{a}_{11}}&{{a}_{21}}&{\cdots}&{ {a}_{p1}}\\{ {a}_{12}}&{{a}_{22}}&{\cdots }&{{a}_{p2}}\\{}&{}&{ \vdots}&{ }\\{ {a}_{1p}}&{{a}_{2p}}&{\cdots}&{ {a}_{pp}}\end{array}\right] $

(9)

1.3.2 KPCA方法原理

KPCA的核心思想就是引入一个非线性函数$ \phi \left(x\right) $，将原始样本集合$ \boldsymbol{X}=\{{\boldsymbol{x}}_{1}, {\boldsymbol{x}}_{2}, \cdots , {\boldsymbol{x}}_{n}\} $从样本空间映射到更高维的特征空间$ \boldsymbol{H} $中^[19]。图 2给出了KPCA方法的基本思路，不同形状代表不同类别的样本。其中，映射函数为

$ \begin{array}{l}\phi :{I}^{m}\to \boldsymbol{H}, {\boldsymbol{x}}_{r}\to \phi \left({\boldsymbol{x}}_{r}\right)\\ r=\mathrm{1, 2}, \cdots , n\end{array} $

(10)

在特征空间$ \boldsymbol{H} $内再使用线性的PCA对样本集合$ \phi \left({\boldsymbol{x}}_{1}\right), \phi \left({\boldsymbol{x}}_{2}\right), \cdots , \phi \left({\boldsymbol{x}}_{n}\right) $进行数据降维，这样问题就变成在特征空间内求解样本协方差矩阵的特征值与特征函数。

特征空间内，样本的协方差矩阵可以表示为

$ {C}^{\mathrm{H}}=\frac{1}{n}\sum\limits_{r=1}^{n}\phi \left({\boldsymbol{x}}_{r}\right){\phi }^{\mathrm{T}}({\boldsymbol{x}}_{r}{)}^{} $

(11)

用$ \lambda $表示协方差矩阵$ {C}^{\mathrm{H}} $的特征值，$ \boldsymbol{v} $表示对应的特征向量，则有

$ {\boldsymbol{C}}^{\mathrm{H}}\boldsymbol{v}=\lambda \boldsymbol{v} $

(12)

定义核函数

$ {\boldsymbol{K}}_{s, j}=〈\phi \left({\boldsymbol{x}}_{s}\right)\times \phi \left({\boldsymbol{x}}_{j}\right)〉 $

(13)

特征向量v在$ \phi \left({\boldsymbol{x}}_{1}\right), \phi \left({\boldsymbol{x}}_{2}\right), \cdots , \phi \left({\boldsymbol{x}}_{n}\right) $构成的空间内，因此必定存在系数集$ {a}_{s}(s=1, \cdots , k) $，使得

$ \boldsymbol{v}=\sum\limits_{s=1}^{k}{a}_{s}\phi \left({\boldsymbol{x}}_{s}\right) $

(14)

式中k表示组成系数集的子集的数量。

将式(13)和式(14)代入式(12)中，得到

$ k\lambda \boldsymbol{K}\boldsymbol{a}={\boldsymbol{K}}^{2}\boldsymbol{a} $

(15)

式中$ \boldsymbol{a}=[{a}_{1}, \cdots , {a}_{k}{]}^{\mathrm{T}} $。因此式(15)等价于

$ k\lambda \boldsymbol{a}=\boldsymbol{K}\boldsymbol{a} $

(16)

需要指出的是，虽然假设原始数据已经进行过中心化处理，但经过映射后的数据并非零均值化的结果，因此用$ \tilde{\boldsymbol{K}} $代替$ \boldsymbol{K} $，对其进行中心化，得到

$ \tilde{\boldsymbol{K}}=\boldsymbol{K}-{\boldsymbol{l}}_{n}\boldsymbol{K}-\boldsymbol{K}{\boldsymbol{l}}_{n}+{\boldsymbol{l}}_{n}\boldsymbol{K}{\boldsymbol{l}}_{n} $

(17)

式中$ {\boldsymbol{l}}_{n} $为$ n\times n $的矩阵且其中元素值都为$ 1/n $。

1.4 模型试算

使用模型对子波整形方法进行测试，模型选取主频为50 Hz的Morlet子波，输入后进行子波整形，与同极性的反射系数序列合成为地震记录。处理前、后的合成记录如图 3所示，显然，整形后的子波具有更窄的主瓣，调谐现象变弱，能更好地与真实反射系数序列对应，具有较高的分辨率。

为进一步验证方法对薄层的识别能力和降噪效果，利用表 1参数构建楔状地质模型。其中，地层②为楔状层。使用25 Hz的Ricker子波作为震源，得到原始地震记录，添加45 dB的高斯噪声后，得到如图 4a的模型记录。

对模型数据进行匹配追踪处理，选出20、25、30 Hz的剖面数据，进行子波整形，整形后的数据经KPCA方法处理，得到楔状模型记录(图 4b)。

模型中实际尖灭点位于第22道，合成记录中调谐作用出现在第49道(图 4a中红色虚线位置)，并将其作为识别出的尖灭点位置；由图 4b可见，尖灭点位置出现在第42道。同时，处理后数据增强了弱反射，并剔除了大部分高斯噪声，有利于识别薄层。

2 实际应用

现有M区实际模型数据，该区储层整体较薄，平均厚度为17.2 m，主要为颗粒泥晶灰岩和颗粒灰岩，具有多层叠置、横向变化快的特点。同时，该区断层非常发育，有12条主要断层，断层产状较陡。常规的地震勘探技术不能满足油藏描述要求，增加了勘探开发的风险。

因此，本文以该区地震资料为基础，采用基于匹配追踪与KPCA的高分辨率处理方法，从断层识别和薄层刻画的角度，对地震资料开展智能重建研究，以期提高资料分辨率。

2.1 断层识别

在地震勘探过程中，细微断层会导致注水效果变差，影响采油效果^[20]。但受地下复杂构造以及地震资料分辨率的影响，很难分辨地层中的细微断层及微小的错动^[21-22]。在图 5a原始地震剖面中，红色椭圆部分解释为同相轴变形，但细微的错动很难判别是否为断层；经过重建后(图 5b)，可以看出相对于原始数据的连续结构，地层存在明显的位移，同相轴发生错断，可以证实为细微断层，蓝色线条为断层解释结果。

2.2 薄层刻画

薄互层储集层是油气藏勘探的重要目标，当层的厚度较小、间隔较近时，同相轴受到相邻地层的影响而被掩盖^[23]。在资料解释过程中，很难分辨薄互层的底部和顶部界面反射，两个薄层的同相轴相互干扰，从而呈现低频、强振幅反射特征，很容易导致错误的解释。

在原始地震剖面（图 6a）中，T10层位的“粗壮”同相轴之间十分模糊，聚焦性较差，薄层分辨率很低；经过本文方法处理后，刻画的层位特征更准确。在重建剖面（图 6b）的紫色区域，原本难以判别的薄层变得清晰可辨，同相轴连续完整，伽马测井曲线呈低值响应特征，有效揭示了厚度约为16 m的砂岩薄层顶底界面。因此，高分辨处理后的剖面具有较强的井震相关性且数据可靠。

2.3 频谱分析

在地震记录中，地震信号的能量会聚焦在一定的时频范围内，而噪声信号往往分布在整个时频空间中^[24]。在重建信号的过程中，选择部分优势频段，提取原始地震信号的有效频段，使信号能量集中在优势频段内，更能凸显目标储层的构造特征。图 7b为高分辨率处理后的井震标定剖面频谱分析结果，相较于原始信号(图 7a)，信号能量集中在优势频带附近，分辨地层的能力得到明显提升，有利于特殊地质构造的识别和储层的精准预测。

3 结论

本文将匹配追踪算法和核主成分分析方法应用于地震信号重建，采用子波整形技术对匹配追踪后的子波进行旁瓣抑制，在保证信号保真度的情况下压制了噪声。通过核主成分分析方法重建地震信号，最大程度凸显了信号的主要特征，信号叠加的同相轴更清晰，分辨率得到明显提升。

实际资料应用表明，该方法适用于地震信号的高分辨率处理。经过重建信号，对目标储层的分辨能力得到提升，对地层的构造特征刻画更精细，能有效识别细微断层、薄互层等较难判别的地质现象，为后续的地震资料解释、储层预测提供了可靠的依据。

在实际应用过程中发现方法仍存在一定的局限性，主要体现在以下几个方面：①本文的研究目的主要是针对细微断层和薄层的细节刻画，并未验证对其他地质构造的处理效果；②匹配追踪及子波整形对较大断层响应较好，经过KPCA重建后对小断层的响应更敏感；③该方法虽然提高了薄层的分辨率，但重构信号时选取部分优势频段会在一定程度上丢失一些细节，针对资料保真、保幅仍需进一步攻关研究。