1. 引言

地震物理勘探法是地球物理学中的一种获取如油气、矿产等能源的地质结构具体分布信息的常用方法。地震波在地层传播过程中不可避免会受到大地低通滤波的影响，因而造成高频信息严重衰减，使得地震信号中的有效信息大大减少。近几年来，提高地震资料的分辨率、增加地震信号的信噪比逐渐成为研究热点[1]。反褶积算法是提高地震资料分辨率的重要方法之一，为了获取更加准确的地震资料，各种各样的反褶积方法不断被提出来。反褶积方法是通过压缩地震记录中的地震子波，压制交混回响与多次波，从而提取反射系数的一种数据处理方法[2]。传统的反褶积方法有，维纳滤波(Wiener)滤波、f-k滤波、最小平方反褶积、预测反褶积、同态反褶积、地表一致性反褶积与反Q滤波，以及本文会涉及的脉冲反褶积与L1范数稀疏脉冲反褶积[3]。近年来，一些新的方法被提出来对反褶积方法进行改进和完善，例如基于压缩感知理论的反褶积方法，自适应的迭代阈值收缩法、由数据驱动的深度褶积网络模型等。

1942年，维纳(N. Wenier)最先提出维纳滤波，维纳滤波就成为了地震数据处理中广泛使用的数字滤波方法之一[4]。同时，根据最小平方原理，改变滤波器的期望输出，可以推出多种反褶积方法，如最小平方反褶积、脉冲反褶积等。1954年，Robinson研究表明地震信号可以被预测分解，从而提出了预测反褶积方法。1967年Robinson讨论了预测反褶积和维纳滤波之间的关系，认为预测反褶积实质上也属于最小平方反褶积[5]。Peacock和Treitel (1969)把预测反褶积应用于实际地震资料处理，结果表明该方法能够有效压制多次波和鸣震[6]。Ulrych (1971)将同态系统引入反褶积方法，提出了同态反褶积[7]。他认为同态反褶积不需要假设地层反射系数在统计上是白噪、地震子波是最小相位子波，其在复赛谱域使用滤波技术就能把地震子波与反射系数分离。近年来，基于深度卷积神经网络和基于模态分解的反褶积方法不断被提出。张联海(2021)等提出由数据驱动的深度卷积神经网络DCNN模型无需再次设置正则化参数即可用于求解稀疏反褶积问题[8]。所提DCNN模型还采用多分辨率分解和残差学习等技术以提高网络的表达能力。倪文军(2023)等设计了子波与反射系数串联、迭代、更新的实现方案，利用已知模型进行网络的预训练，将基于目标数据体提取的有效子波作为靶区数据反褶积的初始子波，进行子波整形反褶积处理[9]。邵佳(2023)等结合EMD反褶积技术与构造导向约束滤波技术提高三维地震资料的分辨率，利用构造导向约束滤波技术滤除EMD反褶积资料的部分噪声[10]。卫泽(2023)等针对同态域子波与反射系数不易彻底分离的问题，提出在同态域使用变分模态分解(Variational Mode Decomposition, VMD)方法分离子波与反射系数信息，提高了同态反褶积的抗噪能力[11]。

基于维纳滤波的反褶积方法需要在一定的假设条件下才能取得相对理想的结果。最终的反褶积处理结果与滤波因子直接相关，原始数据在处理过程中并没有起到直接的约束作用。预测反褶积广泛用于压制多次波，但其预测步长等参数的选择很困难。脉冲反褶积通常有两个前提：一是地震子波为最小相位子波，二是反射系数序列是白噪序列。脉冲反褶积的处理结果容易受到噪声的干扰，应用时需要多加注意。L1范数，L2范数正则化约束的反褶积方法中，滤波因子的作用至关重要，微小的噪声干扰都会引起处理结果的巨大变化[12]。

本文首先介绍了反褶积、维纳滤波、最小平方反褶积、脉冲反褶积和稀疏脉冲反褶积的基本概念。并通过模拟了反射系数序列，利用最小相位子波合成地震记录，最后用脉冲反褶积与稀疏脉冲反褶积作用于地震数据来观察反褶积方法的作用效果。传统的反褶积方法都要以压缩子波为核心，然而稀疏脉冲反褶积将反褶积方法带向了一个新方向，也就是观测反射系数幅值以及位置信息，从而突破了有效信号频带的限制，极大提高地震信号的分辨率。最后利用数据，考查两种反褶积方法的作用效果，挑选出某些道出来观察作用效果，利用均方误差来考察反褶积的效果[13]。

2. 脉冲反褶积的方法原理

反褶积是地震资料提高分辨率和信噪比最常用的处理方法之一，叠前叠后都可以使用。反褶积的主要作用是压缩地震子波，提高地震子波一致性，进而提高地震资料的分辨率和信噪比，减少噪音，从而提高地震资料解释的精度[14]。

2.1. 反褶积

2.1.1. 地震记录的褶积模型

1、理想模型

若地震波以脉冲 $b\delta \left(t\right)$ 形式激发，经过地层时无吸收、透射和多次反射等因素的影响，而且整个传播过程不存在随机干扰，检波器记录下地表质点随时间的振动图像，即地震记录或者地震道

$x\left(t\right)=b\left(t\right)\ast \xi \left(t\right)$ (1)

这里*表示褶积运算， $x\left(t\right)$ 为地震记录， $\xi \left(t\right)$ 为反射系数序列， $b\left(t\right)$ 指的是地震子波。这时得到的实际上是反射系数序列，同实际的地震记录相比，具有更高的分辨率，高频信息更加丰富[2]。

2、实际模型

实际地震记录由有效波 $s\left(t\right)$ 和干扰波 $n\left(t\right)$ 组成。且有

$x\left(t\right)=s\left(t\right)+n\left(t\right)$ (2)

这里的有效波主要指一次反射波，一次反射波可由以下褶积模型表示。

$s\left(t\right)=b\left(t\right)\ast \xi \left(t\right)$ (3)

2.1.2. 地震记录的反褶积模型

反褶积主要通过数学方法，提高地震记录分辨率，从而更加接近反射系数剖面，得到地下介质精确的物理性质与结构。在假设地震记录不含干扰的情况下，由(2)、(3)两式可得：

褶积过程： $x\left(t\right)=b\left(t\right)\ast \xi \left(t\right)$ 。

对应频率域形式滤波过程： $X\left(\omega \right)=b\left(\omega \right)\Xi \left(\omega \right)$ 。

令

$A\left(\omega \right)=\frac{1}{B\left(\omega \right)}$

则可得到反滤波过程： $\Xi \left(\omega \right)=A\left(\omega \right)X\left(\omega \right)$ 。

写成时间域形式反褶积过程： $\xi \left(t\right)= a\left(t\right)\ast x\left(t\right)$ 。

2.1.3. 反褶积问题的特点

1、反褶积结果存在多解性

对于反褶积方程 $\xi \left(t\right)= a\left(t\right)\ast x\left(t\right)$ ，一般情况下只有地震记录 $x\left(t\right)$ 已知的。而另外两个函数是未知的，因此存在多样解。

2、分辨率与信噪比相互制约

反褶积的目的是得到反射系数(频谱很宽，接近白谱)。但是实际的地震资料不可能没有噪音，反褶积把分辨率提高的同时，把有些频段(高频段和低频段)的噪声也放大了，使信噪比下降。采用均方根误差(Root Mean Square Error, RMSM)来衡量反褶积效果[15]

$\text{RMSE}=\sqrt{\frac{1}{nt\times nx}\underset{i=1}{\overset{nt}{{\displaystyle \sum }}}\underset{j=1}{\overset{nx}{{\displaystyle \sum }}}{\left({\xi }_{i,j}-{\hat{\xi }}_{i,j}\right)}^2}$ (4)

式中， ${\hat{\xi }}_{i,j}$ 代表第i个采样点、第j道的真实反射系数，nt和nx分别代表采样点数和地震导数。RMSE越小，意味着反褶积效果越好[16]。

2.2. 最佳维纳滤波

最佳维纳滤波通过滤波器实际输入与期望输出的误差平方和最小从而计算滤波因子。基本原理如下：

1、求解关系

输入信号： $x\left(t\right)=\left(x\left(0\right),x\left(1\right),\cdots ,x\left(n\right)\right)$

滤波因子： $h\left(t\right)=\left(h\left(0\right),h\left(1\right),\cdots ,h\left(m\right)\right)$

实际输出：

$y\left(t\right)=h\left(t\right)\ast x\left(t\right)=\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\ast x\left(t-\tau \right)=\left(y\left(0\right),h\left(1\right),\cdots ,h\left(M\right)\right),M=m+n$

期望输出： $d\left(t\right)=\left(d\left(0\right),d\left(1\right),\cdots ,d\left(m\right),d\left(M\right)\right)$

输出误差： $e\left(t\right)=d\left(t\right)-y\left(t\right)$

误差能量：

$Q\left(t\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}{e}^2\left(t\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}{\left[d\left(t\right)-y\left(t\right)\right]}^2=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}{\left[\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\ast x\left(t-\tau \right)-d\left(t\right)\right]}^2$

为了求取合适的滤波因子 $h\left(t\right)$ ，采用最小平方法，即令误差平方和(总能量误差) Q为最小。这样，滤波因子的求解问题就转换为一个误差能量Q对滤波因子 $h\left(t\right)$ 求极值的问题。

将Q对滤波因子 $h\left(t\right)$ 求偏导，并令其等于零

$\begin{array}{c}\frac{\partial Q}{\partial h\left(s\right)}=\frac{\partial }{\partial h\left(s\right)}\underset{t=0}{\overset{M}{{\displaystyle \sum }}}{\left[\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\ast x\left(t-\tau \right)-d\left(t\right)\right]}^2\\ =\underset{t=0}{\overset{M}{{\displaystyle \sum }}}\frac{\partial }{\partial h\left(s\right)}{\left[\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\ast x\left(t-\tau \right)-d\left(t\right)\right]}^2\\ =2\underset{t=0}{\overset{M}{{\displaystyle \sum }}}\left[\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\ast x\left(t-\tau \right)-d\left(t\right)\right]x\left(t-s\right)=0\end{array}$

由此得出

$\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}h\left(\tau \right)\underset{t=0}{\overset{M}{{\displaystyle \sum }}}x\left(t-\tau \right)x\left(t-s\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}d\left(t\right)x\left(t-s\right)\left(s=0,1,\cdots ,m\right)$

令

$r_{xx}\left(\tau -s\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}x\left(t-\tau \right)x\left(t-s\right)$

$r_{dx}\left(s\right)=\underset{\tau =0}{\overset{M}{{\displaystyle \sum }}}d\left(t\right)x\left(t-s\right)$

即有

$\underset{\tau =0}{\overset{m}{{\displaystyle \sum }}}{r}_{xx}\left(\tau -s\right)x\left(t-s\right)=r_{dx}\left(s\right)$ (5)

上式 $r_{xx}\left(\tau -s\right)$ 是延迟为 $\tau -s$ 的输入信号自相关， $r_{dx}\left(s\right)$ 是延迟为s的输入信号与期望输出的互相关。则上式可写成矩阵形式

$\left[\begin{array}{cccc}{r}_{xx}\left(0\right)& r_{xx}\left(1\right)& \cdots & r_{xx}\left(m\right)\\ r_{xx}\left(1\right)& r_{xx}\left(0\right)& \cdots & r_{xx}\left(m-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ r_{xx}\left(m\right)& r_{xx}\left(m-1\right)& \cdots & r_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}a\left(0\right)\\ a\left(1\right)\\ \begin{array}{c}⋮\\ a\left(m\right)\end{array}\end{array}\right]=\left[\begin{array}{c}{r}_{dx}\left(0\right)\\ r_{dx}\left(1\right)\\ \begin{array}{c}⋮\\ r_{dx}\left(m\right)\end{array}\end{array}\right]$ (6)

通过(6)式，可以得到滤波因子 $a\left(t\right)$ ，再与输入信号 $x\left(t\right)$ 褶积，最后可得输出信号 $y\left(t\right)$ 上式中， $x\left(t\right)$ 自相关矩阵为Toeplitz矩阵。

2.3. 最小平方反褶积

将最佳维纳滤波原理应用于反褶积，即为最小平方反褶积。地震子波 $b\left(t\right)$ 相当于输入信号 $x\left(t\right)$ ，期待输出 $d\left(t\right)$ 可以是任意期望的波形。 $b\left(t\right)$ 未知，无法得知Toeplitz矩阵，下面推导子波自相关与地震道自相关的关系[16]。

褶积模型可以表示为： $x\left(t\right)=b\left(t\right)\ast \xi \left(t\right)$ ；

上式在频率域可以表示为： $X\left(\omega \right)=B\left(\omega \right)\Xi \left(\omega \right)$ ；

即有

${\left|X\left(\omega \right)\right|}^2={\left|B\left(\omega \right)\right|}^2{\left|\Xi \left(\omega \right)\right|}^2$

如果假设反射系数为白噪，其自相关函数为：

${\varphi }_{rr}\left(t\right)=\{\begin{array}{l}e\left(t=0\right)\\ 0\left(t\ne 0\right)\end{array}$ (t为常数)

则有 ${\left|X\left(\omega \right)\right|}^2=e{\left|B\left(\omega \right)\right|}^2$ 。

上式表明：地震道自相关函数与地震子波自相关系数相差一个比列系数e。e为常数，相当于在反滤波因子上乘以一个常数，因此可以用地震道自相关代替子波自相关。

反滤波因子 $a\left(t\right)$ 求解方程如下：

$\left[\begin{array}{cccc}{\varphi }_{xx}\left(0\right)& {\varphi }_{xx}\left(1\right)& \cdots & {\varphi }_{xx}\left(m\right)\\ {\varphi }_{xx}\left(1\right)& {\varphi }_{xx}\left(0\right)& \cdots & {\varphi }_{xx}\left(m-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\varphi }_{xx}\left(m\right)& {\varphi }_{xx}\left(m-1\right)& \cdots & {\varphi }_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}a\left(0\right)\\ a\left(1\right)\\ \begin{array}{c}⋮\\ a\left(m\right)\end{array}\end{array}\right]=\left[\begin{array}{c}{\varphi }_{dx}\left(0\right)\\ {\varphi }_{dx}\left(1\right)\\ \begin{array}{c}⋮\\ {\varphi }_{dx}\left(m\right)\end{array}\end{array}\right]$ (7)

2.4. 脉冲反褶积

有一种特殊情况，当 $d\left(t\right)=\delta \left(t\right)$ 时

$\delta \left(t\right)=\{\begin{array}{l}1\left(t=0\right)\\ 0\left(t\ne 0\right)\end{array}$ .

$r_{dx}\left(s\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}d\left(\tau \right)x\left(t-s\right)=\underset{t=0}{\overset{M}{{\displaystyle \sum }}}\delta \left(\tau \right)x\left(t-s\right)=\{\begin{array}{l}b\left(0\right)s=0\\ 0s\ne 0\end{array}$

方程组可以进一步简化

$\left[\begin{array}{cccc}{r}_{xx}\left(0\right)& r_{xx}\left(1\right)& \cdots & r_{xx}\left(m\right)\\ r_{xx}\left(1\right)& r_{xx}\left(0\right)& \cdots & r_{xx}\left(m-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ r_{xx}\left(m\right)& r_{xx}\left(m-1\right)& \cdots & r_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}a\left(0\right)\\ a\left(1\right)\\ \begin{array}{c}⋮\\ a\left(m\right)\end{array}\end{array}\right]=\left[\begin{array}{c}b\left(0\right)\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right]$ (8)

其中， $b\left(0\right)$ 为一个系数，这样求出的反子波来做反褶积处理，输出统一乘上了一个比例系数。因此上式可以化简为：

$\left[\begin{array}{cccc}{r}_{xx}\left(0\right)& r_{xx}\left(1\right)& \cdots & r_{xx}\left(m\right)\\ r_{xx}\left(1\right)& r_{xx}\left(0\right)& \cdots & r_{xx}\left(m-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ r_{xx}\left(m\right)& r_{xx}\left(m-1\right)& \cdots & r_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}a\left(0\right)\\ a\left(1\right)\\ \begin{array}{c}⋮\\ a\left(m\right)\end{array}\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right]$ (9)

因此，只要知道地震子波的自相关，即可求出反子波，在得出反褶积的输出，由于期待输出为 $\delta$ 脉冲，这种反褶积方式被称为脉冲反褶积。上式使用的前提条件是地震子波为最小相位子波，为了计算稳定，需要于柏华处理，给Toeplitz的对角线元素加上一个常数 $\lambda r_{xx}\left(0\right)$ ，可得

$\left[\begin{array}{cccc}\left(1+\lambda \right)r_{xx}\left(0\right)& r_{xx}\left(1\right)& \cdots & r_{xx}\left(m\right)\\ r_{xx}\left(1\right)& \left(1+\lambda \right)r_{xx}\left(0\right)& \cdots & r_{xx}\left(m-1\right)\\ ⋮& ⋮& \ddots & ⋮\\ r_{xx}\left(m\right)& r_{xx}\left(m-1\right)& \cdots & \left(1+\lambda \right)r_{xx}\left(0\right)\end{array}\right]\left[\begin{array}{c}a\left(0\right)\\ a\left(1\right)\\ \begin{array}{c}⋮\\ a\left(m\right)\end{array}\end{array}\right]=\left[\begin{array}{c}1\\ 0\\ \begin{array}{c}⋮\\ 0\end{array}\end{array}\right]$ (10)

以上方程组可以看出， $\lambda$ 太小，对方程求取稳定解帮助不大， $\lambda$ 太大，反褶积的作用变小。特别地，当 $\lambda \to \infty$ ， $a\left(t\right)\to \delta \left(t\right)/c$ ，其中 $c=\left(1+\lambda \right)r_{xx}\left(0\right)$ 。这时反褶积输出就等于地震记录本身，反褶积根本不起作用。在实际操作中，白噪系数一般取0.5%~5%，最多不超过10%。

脉冲反褶积程序实现过程中应注意主要参数的选择：

(1) 反滤波算子长度的选择：

反滤波算子 $a\left(t\right)$ 的长度m理论上可以任意选择。在一个剖面上，可以通过试验来选择最佳的长度，但是通常情况下选择为与地震子波长度接近，一般为80~200 ms。同时，计算时根据地震数据的采样率换算为采样点数。

(2) 相关时窗长度的选择；时窗长度 $m+n$ 至少应大于反滤波算子长度m的两倍。

综上所述，实现脉冲反褶积主要有以下几个步骤：

选取合适的时间窗和地震道范围的地震数据，对每一道地震数据进行如图1步骤的处理：

Figure 1. Flowchart of the spike deconvolution flow chart

图1. 脉冲反褶积流程图

2.5. 稀疏脉冲反褶积

脉冲反褶积需要假设子波为最小相位，并且反射系数为白噪。如果不满足条件，反褶积效果会不理想，进行脉冲反褶积过程中，提高分辨率的同时会降低信噪比。稀疏脉冲反褶积对地震子波的相位没有要求，只需要加入反射系数序列的先验信息，将反褶积问题转化为极值问题，最后利用数值迭代算法就可以求解[16]。

2.5.1. 稀疏脉冲反褶积原理

首先考虑褶积模型：

$x=b\ast \xi +n$

其中，地震记录表示为x，地震子波表示为b，反射系数表示为 $\xi$ ，n为干扰。矩阵表示为

$X=B\ast \xi +n$

$X$ 为地震列向量， $\xi$ 为反射系数列向量， $B$ 为子波褶积矩阵， $n$ 为随机噪声列向量。若设地震道长度为k，则：

$\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_{k-1}\\ x_k\end{array}\right]=\left[\begin{array}{ccccc}{b}_0& b_{-1}& \cdots & \cdots & b_{-\left(k-1\right)}\\ b_1& b_0& b_{-1}& \ddots & ⋮\\ ⋮& b_1& b_0& b_{-1}& ⋮\\ ⋮& \ddots & b_1& b_0& b_{-1}\\ b_{k-1}& \cdots & \cdots & b_1& b_0\end{array}\right]\left[\begin{array}{c}{\xi }_1\\ {\xi }_2\\ ⋮\\ {\xi }_{k-1}\\ {\xi }_k\end{array}\right]+\left[\begin{array}{c}{n}_1\\ n_2\\ ⋮\\ n_{k-1}\\ n_k\end{array}\right]$

在矩阵 $B$ 已知的情况下，希望找到反射系数 $\xi$ ，即寻找使得下列函数J取最小值时的 $\xi$ 。

$J=\underset{\xi }{\min }{\Vert X-B\xi \Vert }_2^2$ (11)

对上式进行求导，并令导数为0，可得：

$\frac{\partial J}{\partial {\xi }_i}=0$

上式可变形化简为：

$B^{\text{T}}\left(B^{\text{T}}\xi -X\right)=0$

即：

$\xi ={\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}X$ (12)

$B^{\text{T}}B$ 称为协方差矩阵， ${\left(B^{\text{T}}B\right)}^{-1}{B}^{\text{T}}$ 称为最小平方逆。地球物理问题中， $B^{\text{T}}B$ 一般没有逆，上式为一个病态问题，不能得到适定的解。所以需要对模型 $\xi$ 进行正则化约束。目标函数可以改为：

$J=\underset{\xi }{\min }\frac{1}{2}{\Vert X-B\xi \Vert }_2^2+\lambda \phi \left(\xi \right)$ (13)

上式中的 $\phi \left(\xi \right)$ 为正则化项，对模型 $\xi$ 的解做先验条件约束， $\lambda$ 为权衡因子，用来调节目标函数中拟合误差项与正则化项的比重。正则化的选取取决于反射系数满足的先验条件。一般认为反射系数序列由较大的稀疏反射系数和服从高斯分布的微小反射系数背景组成。忽略微小的高斯背景，可认为反射系数具有稀疏性[16]。

2.5.2. 正则项选取及实现

对于(10)式中求解 $X$ 存在不稳定的情况下，因此考虑加入正则项。正则项包括L₁。正则项和L₂正则项，L₁正则项可使某些模型参数为0，用于特征提取；L₂正则项可防止过度拟合。本文的稀疏性采用L₁范数衡量[17]。

因此目标函数的正则项可以设置为：

$\phi \left(\xi \right)={\Vert \xi \Vert }_1$

目标函数可修改为：

$J=\underset{\xi }{\min }\frac{1}{2}{\Vert X-B\xi \Vert }_2^2+\lambda {\Vert \xi \Vert }_1$ (14)

式中， ${\Vert \cdot \Vert }_1$ 表示一范数。

利用匹配追踪算法对加入L₁范数的稀疏脉冲反褶积方法进行求解。匹配追踪算法是一种稀疏分解算法，它的基本原理是通过不断的迭代，从原子库中匹配出最好的原子，然后通过线性叠加组合实现对信号的不断逼近。给定一个信号长度为N的待分解信号 $f\in H$ ，用匹配追踪算法对信号进行稀疏表示的步骤如下：

① 构造恰当的过完备原子库 $D={\left\{g_{\gamma }\right\}}_{\gamma \in \Gamma }$ ：选原子→参数化→离散化→归一化。

利用非零相位Ricker子波建立完备原子库，公式为：

$r^{\prime }\left(t\right)=r\left(t\right)\cos \omega ˚-r^{*}\left(t\right)\sin \omega ˚$

在这里， $r\left(t\right)$ 为非零相位Ricker子波， $\omega ˚$ 为相位旋转角度， $r^{*}\left(t\right)$ 为 $r\left(t\right)$ 的Hilbert变换。上式可写为：

$g_{\gamma \in \Gamma }=r\left(t\right)\cos \omega ˚-r^{*}\left(t\right)\sin \omega ˚$

将非零相位Ricker子波公式离散化， $\Gamma$ 是时频参数的集合。 $f_p$ 为主频参数， $\omega$ 和u为相位参数和位移参数。时频参数按下式离散化：

$\gamma =\left(i\ast d{f}_p,j\ast d\omega ,n\ast du\right)$

对上式进行归一化，令 $\Vert g_{\gamma }\left(t\right)\Vert =1$ 。 $f_p$ 、 $\omega$ 和u确定子波的形状和位置。

② 在过完备原子库中寻找最佳匹配原子 $g_{{\gamma }_0}$ ， $g_{{\gamma }_0}$ 满足：

$\left|\langle f,g_{{\gamma }_0}\rangle \right|=\underset{\gamma \in \Gamma }{\sup }\left|\langle f,g_{\gamma }\rangle \right|$

其中， $\left|\langle f,g_{{\gamma }_0}\rangle \right|$ 是信号f与原子 $g_{{\gamma }_0}$ 的内积绝对值。此时，信号f可表示为：

$f=\langle f,g_{{\gamma }_0}\rangle g_{{\gamma }_0}+{\Vert R^{1}f\Vert }^2$

残差 $R^{1}f$ 越小， $\langle f,g_{{\gamma }_0}\rangle g_{{\gamma }_0}$ 与f越接近，即此时 $\langle f,g_{{\gamma }_0}\rangle g_{{\gamma }_0}$ 为原信号最佳逼近， $g_{{\gamma }_0}$ 为最佳匹配原子。

③ 分解残差信号，不断迭代，直至残差达到一定阈值或者迭代次数达到一定。

第 $k+1$ 次迭代表达式可以表示为： $R^{k}f=\langle R^{k}f,g_{{\gamma }_k}\rangle g_{{\gamma }_k}+R^{k+1}f$ 。

如图2展示了匹配追踪算法对残差进行连续正交分解的过程，当达到一定迭代次数，或者残差达到一定阈值，停止残差信号分解[17]。

Figure 2. The schematic diagram of the $k+1$ iteration of matching pursuit

图2. 匹配追踪第 $k+1$ 次迭代示意图

④ 重构原始信号。

可得，原始信号可以表示为： $f\approx {\displaystyle {\sum }_{k=1}^M\langle R^{k}f,g_{{\gamma }_k}\rangle g_{{\gamma }_k}}$

3. 方法的处理效果

为了测试脉冲反褶积的处理效果，首先合成模拟地震记录。由(1)式可知，地震记录可由反射系数序列与地震子波褶积得到。反射系数序列代表地下不同地层界面介质的反射性质。为了模拟地震记录，先随机生成一个反射系数序列，再让该序列与选择的地震子波进行褶积。利用正态分布随机合成采样间隔为1 ms、长度为900个采样点的反射系数序列。

褶积(Convolution)是将反射系数序列与子波结合以模拟地震记录的过程。通过褶积，可以模拟地震波在地下传播时遇到不同地层界面的反射情况。同时，选择最小相位子波Ricker子波与其褶积合成地震记录。通过褶积模型，利用最小相位子波合成地震记录。分别采用脉冲反褶积方法和稀疏脉冲反褶积进行处理，并且计算处理之后的均方误差。

3.1. 脉冲反褶积方法处理效果

采用脉冲反褶积处理该地震记录，稳定常数设置为0.02，滤波器的长度设置为60，处理结果如图3。

由图4可知，反射系数的强振幅在脉冲反褶积过程中得到较好恢复，但是反褶积过程在提高分辨率的同时，在反褶积结果中包含了不可忽略的随机噪声，反射系数不能完全恢复原样。通过均方误差的计算，脉冲反褶积处理前的均方误差为0.048，处理之后的均方误差下降到0.012。所以，脉冲反褶积可以拓宽地震记录的频带，提高地震数据的分辨率。

Figure 3. The effect spike deconvolution processing

图3. 脉冲反褶积处理效果

Figure 4. The analysis of spiking deconvolution spectrum seismic records

图4. 地震记录脉冲反褶积频谱分析

3.2. 脉冲反褶积方法处理效果

为了测试稀疏脉冲反褶积的处理效果，利用以上合成地震数据。通过褶积模型，与最小相位子波合成地震记录。最后采用L1范数约束的稀疏脉冲反褶积方法对合成地震记录进行处理。处理结果如图5所示。

Figure 5. The effect of sparse spike deconvolution

图5. 稀疏脉冲反褶积效果

由图6可知，该情况下的频谱对比展示于下图。反射系数的振幅在稀疏脉冲反褶积过程中得到较好恢复。但因为噪音不可忽略，一些振幅的反射系数没有很好地恢复过来。通过均方误差的计算，稀疏脉冲反褶积处理前的均方误差为0.086，处理之后的均方误差下降到0.014。所以，稀疏脉冲反褶积可以有效提高地震数据的分辨率。

Figure 6. Analysis of sparse pulse deconvolution spectrum of seismic records

图6. 地震记录稀疏脉冲反褶积频谱分析

4. 结论

地震物理勘探中，由于大地滤波的存在，地震波能量会损耗、衰减，再加之地层介质和人为活动产生的噪音无法避免。提高地震资料的分辨率，对于地震解释具有重要的意义。本文主要以地震数据处理为目标，以反褶积原理作为理论指导，对比了脉冲褶积和稀疏脉冲反褶积对于地震数据处理的效果。进一步了解反射系数与地震数据之间有何映射关系。通过合成模拟数据证明，脉冲反褶积与稀疏脉冲反褶积在地震数据的处理上有一定的作用效果，可以有效地提高地震数据分辨率。

本文虽然在地震数据反褶积中验证了脉冲反褶积和稀疏脉冲反褶积的确在反褶积过程中具有一定效果。但是还缺乏一定的不足，研究方法和思路依然需要不断完善和创新。主要体现在对于脉冲反褶积和稀疏脉冲反褶积中的滤波因子、滤波系数以及其他一些参数的优化选择还不够。同时，关于脉冲反褶积和稀疏脉冲反褶积的求解方法也需要一定的精进。

参考文献