1 引言

勘探开发过程中，可以通过构造解释、属性分析、岩性反演等手段识别储层有利区，但是地震分辨率问题直接关系方法的预测精度.传统上，薄层分辨是分辨率问题的重要任务，常以子波能量展布和震荡长度来定量分析，决定于地震记录的主频和频宽范围(Sheriff，1977；Chung and Lawton, 1995；Schoenberger，1974).大多数高分辨率处理方法也相应地关注于分辨薄层和提高频宽，但都面临着信噪比与分辨率之间的矛盾(李庆忠，1993；Robinson and Treitel, 1980)，或处理结果不稳定(黄储德，1992)，或提升效果不明显(Lu and Liu, 2007)，单纯讨论方法的信噪比和分辨率是片面的.现在，越来越复杂的勘探对象，比如岩性油气藏和非常规油气的勘探和开发，要求高分辨率算法不仅仅分辨薄层，更需要提供更丰富的波形和相位信息(Lu，2009).

理论表明，反射系数的频谱形式为有限采样长度的正弦叠加，基于此可以将对反射系数或者薄层厚度估计问题转为对叠加正弦的谱分析(Marfurt and Kirlin, 2001).Puryear和Castagna(2008)建立了利用地震道的频谱信息计算地层厚度的方法，并表明其能够适用于小于调谐厚度的薄层；同样，基于稀疏先验信息的反演类方法也被应用到频率域，用于高分辨地恢复反射信息(Ulryth，2001)；考虑到反问题的不适定性以及某些脉冲类反演会丢失相位信息(Yuan and Wang, 2013；Velis，2008)，源于非参数谱估计理论的方法开始被应用到高分辨率处理.其中自适应滤波方法，简称为APES(Stoica and Moses, 1997)，主要对叠加正弦信号的幅度和初始相位进行估计：对每个时间位置，算法设计相应的自适应滤波器压制其他位置反射和噪声带来的干扰，以达到对该位置反射幅度及相位的稳定估计；实践表明，自适应滤波方法在提高地震资料分辨率的同时，处理结果相应表现出更丰富的信息(Guo and Wang, 2012).

作为新的地震资料自适应高分辨率处理算法，方法特性及适用条件还缺乏深入讨论，特别针对信噪比和分辨率的关系问题还需要进一步分析.因此，本文在讨论自适应滤波方法理论特性的基础上，根据地震数据模型的特点提出了加权改进算法，并特别对方法改进前后在分辨率和信噪比方面的应用效果进行定量的刻画，进一步讨论了算法的薄层分辨能力.最后，应用典型构造模型的波动方程正演模拟对算法进行验证，并对某探区的三维实际资料进行处理实践.

2 方法原理 2.1 地震采样模型

由于子波带限，地震记录中高频成分以及某些低频成分被压制，出于稳定性和信噪比考虑，基于褶积模型的线性反滤波方法只能恢复有效频带里的信息(Kay and Marple, 1981).从这种意义上说，重建反射系数序列的问题是“病态”的，对垂直分辨率的改善也收效甚微.考虑对地震信号模型重新定义，比如将反射信息参数化为时间位置和反射幅度，可以降低实际需要恢复的信息量，使问题变清晰.自适应滤波方法(APES)则以地震信号的采样模型为基础来处理问题：设反射系数r(t)为不同位置t_k和幅度a_k的K个δ脉冲组成的脉冲串，经地震子波w(t)滤波褶积后为连续地震信号x(t)，然后进行间隔为Δt的离散采样得到记录X(nΔt).这一过程在时间域和频率域的形式分别为(Unser，2000)

(1)

(2)

其中，*代表连续褶积，A/D代表离散采样过程，k代表K个反射脉冲的序号.变量t为连续时间，nΔt为t的离散采样，n代表地震记录采样点序号，N为采样点数，m代表泊松原理中频谱混叠的周期叠加序号，f和nΔf分别为连续和离散采样的频率，为子波频谱，为地震记录频谱.

从采样模型(公式(1)、(2))可以看出，子波的滤波作用为地震反射系数的重建提供了一个条件，也摆出了一个限制：一方面由于地震记录的有限频带宽度( < 1/Δt)，因此离散采样后没有超过1/Δt的高频能量与有限带宽内的能量相叠加，即在一个叠加周期内(1/Δt)混叠效应(公式(2))对有效带宽内的频谱没有影响，等于反射系数的谱乘以子波的谱(Unser，2000；Guo and Wang, 2012)，这就给从频率域寻求估计反射系数的方法提供了前提条件；另一方面，子波的滤波作用压制了反射系数在有效带宽外的能量，有限带宽一方面决定了时间域地震信号的分辨率特征，另一方面给有效恢复反射系数的高频信息带来限制(李庆忠，1993).

2.2 自适应滤波方法(APES)

对地震记录和子波分别进行离散傅里叶变换可以得到频谱序列和，再用子波的谱除地震记录的谱就可以得到反射系数的部分谱序列，其形式为有限长度正弦信号的叠加，形式如下：

(3)

其中，e(jΔf)代表噪声，Δf=1/(NΔt)为频率的离散采样间隔，L和U对应地震数据有效频带的上下限.选取任意时间位置t₀，自适应滤波方法在估计位置t₀的反射幅度a(t₀)时，首先将频谱视为两部分之和：正弦a(t₀)e^{－i2πt₀jΔf}(a(t₀)特指与t₀的对应关系)和由其他的反射正弦和噪声构成的相对干扰，即将式(3) 重新写为

(4)

为了估计对应t₀的反射幅度a(t₀)，方法设计相应的自适应滤波器来压制干扰.首先将序列分成J个M维子向量如下(J=U－L－M+2，T代表矩阵实转置)：

(5)

然后设计M维滤波器H=[h₀, h₁, …h_M-1]_M×1^H([]^H代表矩阵复共轭转置)，要求其尽量压制相对干扰v(nΔf)的同时，正弦分量序列a(t₀)e^{－i2πt₀jΔf}无损地通过，于是自适应滤波器的设计为如下目标函数优化问题：

(6)

这里||代表复数求模，α(t₀)=[1, e^-i2πt₀Δf, …e^{-i2πt₀(M-1)Δf}]_M×1^T.

由公式(6) 可得：

(7)

这里.

为使得目标函数取最小值，因此有a(t₀)=H^Hg(t₀)，而令，则优化问题(公式(6))转化为.H^Hα(t₀)=1.这一优化问题可由柯西施瓦茨原理得到(Stoica et al., 1999)

(8)

对应的.从方法原理可以看出，APES方法是对任意选定的任意时间位置t₀计算自适应滤波器来压制噪声和其他位置反射组成的相对干扰，因此应用方法进行处理时，可以计算每个采样时间位置的反射幅度，从而获得连续采样的波形结果，且结果相对于原地震记录表现出更高的分辨率(具体APES的处理流程请参看文献Guo and Wang, 2012).

2.3 加权的自适应滤波方法(WAPES)

仔细分析APES方法的应用过程可知，原有的滤波器设计优化中没有考虑频率域噪声的具体分布特征，特别是谱除子波后获得的反射系数频谱在低频和高频段的噪声能量会显著提高，这样一来，原方法在提高分辨率的同时，实际降低了信噪比.因为在不同频段噪声水平不一致，所以对滤波器的设计优化过程应该将噪声的能量水平考虑进去，建立优化目标函数时按照不同频段的噪声水平选取相应的权重，从而获得更优的滤波器设计结果.因此，本文根据不同频段反射系数部分频谱的信噪比特征，将自适应滤波器的设计优化问题(公式(6))进行加权可得

(9)

进一步地，推导可得

(10)

这样优化问题(公式(9))的解与常规方法有着相同的形式

(11)

而这里的.假设原地震记录的噪声能量在有效频带内均匀分布，则考虑谱除子波对噪声在不同频段的放大效应，权重因子W_j可以根据子波的振幅谱自适应地进行选择.比如对于离散的子波的振幅谱，权值可以定义为(如图 1所示)

(12)

如果噪声的能量在有效频带内并不满足均匀分布，则需要同时考虑噪声的能量分布在谱除后的信噪比的影响.假设谱除子波前有效频带范围内噪声能量分布的离散形式为{P_n(jΔf)}，则权值可以定义为

(13)

实际上，若考虑频率域信噪比的定义(Zhao Yan, et al., 2014)：

(14)

其中P_s(f)和P_n(f)分别为信号的噪声在频率域的能量分布.假设谱除子波后有效频带内信噪比的离散形式为{PSNR(jΔf)}，则相应的权值定义为

(15)

如果考虑地震记录信噪比或者子波随时间发生变化，则可以采取分时段的方式应用方法进行高分辨率处理：先应用时频域信噪比估算方法(Zhao Yan, et al., 2014)计算每个时段的频率域信噪比，然后计算得到的时变权值应用于分时段的加权处理方法(WAPES).

3 自适应滤波理论实验

针对自适应滤波的理论特性进行分析，这里首先设计简单模型进行实验：如图 2所示，两个反射分别位于0.05 s和0.062 s，幅度一致，与25 Hz主频雷克子波褶积合成记录，可以看出反射对是不可分辨的.反射系数的频谱为对应正弦的叠加.这里选取有效频带为5~65 Hz，间隔1 Hz，分别对位置0.025 s，0.05 s，0.056 s和0.062 s应用方法来考察方法在自适应滤波方面的特性(滤波器长度选择为频带长度的一半).

没有噪声的情形下，如图 3所示，0.025 s处对应的滤波器的幅度效应在两个已知反射位置处为0，在0.025 s处为1，即视两个反射为相对干扰进行压制的同时保持了0.025 s处对应的正弦能量，由于模型在该位置没有反射，所以滤波输出为0，相应的反射幅度估计也为0；对于第一个反射位置0.05 s处，对应的滤波器则将0.062 s的反射视为相对干扰，幅度效应为0，在0.05 s处为1，滤波输出恰为频率对应时间0.05 s的正弦，振幅为反射幅度，从而可以估计出0.05 s处的反射幅度；对于0.062 s，与0.05 s处同理，滤波输出恰为频率对应反射时间，振幅对应反射幅度的正弦；而对于位于两个反射中间的时间点0.056 s，设计的滤波器压制两侧的反射，滤波后结果为0，估计结果在两个反射之间存在波谷凹陷，即薄层反射从视觉上可以被分辨.

考察有噪声的情况，仍然考察同样4个位置对应的自适应滤波器效应.如图 4所示，各时间点滤波器的幅度效应与没有噪声的情形基本相似，但由于噪声的存在，整体的幅度偏低.这可以解释为滤波器的设计过程中除了需要压制其他反射的相对干扰能量之外，还要考虑整体压制噪声的能量.因为噪声是遍布整个记录时间的，所以滤波器幅度效应整体偏低.也就是说，在估计某一时间位置反射幅度的时候，自适应滤波器的设计权衡了噪声压制和其他反射相对干扰压制之间的关系.

需要注意的是有噪声的情况下0.056 s处对应的自适应滤波器对于干扰的压制作用明显变弱，所以滤波输出也为两个未完全压制反射带来的相对干扰和噪声的组合.事实上，随着噪声水平的提高，在估计薄层反射之间的反射幅度之时，自适应滤波器的设计在压制噪声能量和相对干扰之间的权衡效应导致对于薄层干扰的压制能力变弱.由于噪声和干扰压制的不完全，方法估计结果的波形变宽，表现出对薄层的分辨能力变差(图 5).

4 信噪比与分辨率

从对自适应方法理论特征的分析可以看到，由于自适应滤波的权衡效应，噪声会影响自适应滤波方法的分辨率及处理结果的稳定性，因此之后将重点对方法在分辨率和信噪比方面的表现进行讨论，并且对这两方面的应用效果进行定量的刻画.定量分析过程分为三个部分，如图 6所示.

4.1 分辨率定量实验

以某口实测探井为例，利用声波及密度资料计算井旁反射系数，并根据采样模型合成井旁地震记录.在无噪的情形下，自适应方法对35 Hz子波合成记录的处理结果与50 Hz子波的合成记录在分辨细节表现方面基本一致(图 7a)，频谱特征也基本相同(图 7b)，即方法有效地提高了35 Hz记录的主频和频带宽度.为了从频谱特征来对方法的分辨率特征进行定量的刻画，这里借用用于刻画地震资料分辨率特征的两个关键参数：地震资料主频f_m以及半有效带宽f_b，分别定义为(Barnes, 1993a, 1993b)：

(14)

其中P(f)为地震记录的振幅谱，而P(f)²为能量谱.

考虑噪声的影响，选取5至60 Hz主频雷克子波合成记录(步进5 Hz)，分别在0%，1%，5%的噪声情况下应用方法，并通过统计试验来讨论处理结果的主频和带宽：从对主频参数的统计试验结果(图 8a，b)可以看到，改进前后的方法都提高了资料的主频，两者基本呈正相关；而从有效频带宽度统计试验结果(图 8c，d)可以看到方法具有一定程度拓宽频带的能力.而对比原算法(APES)和加权算法(WAPES)在主频，频宽定量效果，可以看到改进方法在拓频能力方面(主频提高，频带拓宽)略低于原方法，这主要源于两个原因：一方面加权方法对高频信息的利用要低于原方法，所以拓频能力略低，这一点从无噪声的统计结果可以看出；另一方面，高频噪声会引起模型道的带宽变化，同时原方法提高了更多的高频噪声能量，所以利用公式估计出的主频和带宽表现为偏高的假象.因此，单纯地从主频和带宽去判断高分辨率方法的处理效果是片面的，还需要从信噪比来考察所拓宽的高频能量的有效程度.这一点在后面对信噪比定量的讨论中会给予详细阐述.

4.2 信噪比定量实验

地震资料的信噪比有很多的衡量方式，这里采用刘洋和李承楚(1997)提到的方法，从频率域来计算信噪比.对一定时窗内的Nx道地震记录在有效频带内的频谱组成的矩阵应用SVD分解得到Nx个奇异值{λ_k}_{k=1, 2, …, Nx}，对应估算的信噪比为

(15)

仍然选取实测井资料进行试验，分别选取25 Hz和33 Hz子波合成10道叠前道集记录(采样间隔1 ms).对两者附加相同的高斯白噪声(5%)，并对25 Hz子波合成的叠前道集分别应用APES和WAPES方法进行处理.如图 9所示，对处理前后的道集分别计算信噪比曲线和频谱并进行对比，发现高频(33 Hz)合成地震道的信噪比要低于低频(25 Hz)，而方法改进前后的处理结果主频都达到了33 Hz左右，但是APES的处理结果的信噪比相对于33 Hz合成记录要低很多，这意味着处理过程中噪声能量被过度放大了，而WAPES则很好地控制了噪声能量，与33 Hz合成记录基本一致.

选取不同主频的地震子波合成叠前道集(10~50 Hz，10 Hz步进)，在所附加噪声完全相同的情形下分别对原始道集，处理结果，及处理结果对应主频的合成叠前道集的信噪比进行估计(关注于0.3 s时间位置).如图 10所示，随着主频的增加，合成记录的频带变宽，而对应的有效频带内的信噪比明显降低.即高分辨率处理方法在提高资料主频、拓宽频带的同时，相对于有效信号，噪声的能量实际上被严重的放大了，相应信噪比也相应的表现为大幅度降低.可以看到改进后算法的信噪比与其处理结果对应主频的合成记录的信噪比相一致，这说明改进的算法拓宽频带的同时噪声能量并没有被过度地放大.

4.3 薄层分辨定量实验

由于主频和带宽参数不能直观刻画方法在薄层分辨方面的特性，且容易受到高频噪声的影响，这里借用Widess(1973)对于薄层反射系数对模型的研究思路，在对方法直观分辨薄层能力的研究中提出了如下的公式：

(16)

其中P(t)表示时间位置t估计所得的振幅值，t₁，t₂正好对应于反射系数对的两个时间位置，P[(t₁+t₂)/2]对应反射系数之间的时间位置所估计的波形幅度.在讨论方法分辨率行为的时候，对于估计结果满足公式(16) 的情形，我们认为薄层反射是可分辨的(如图 11所示).

为了定量讨论方法的薄层分辨能力，我们选取15、25、35 Hz以及45 Hz的雷克子波与不同时间厚度的薄层反射对褶积来合成地震记录(薄层厚度1~40 ms，步进1 ms)，并附加1%的白噪声，分析方法针对不同主频不同薄层厚度的合成记录进行处理的分辨率表现(图 12a)，实验针对每种情形各实现500次统计试验，统计方法在相应主频下对不同时间厚度的薄层反射可分辨的百分比.然后选取30 Hz的雷克子波与不同时间厚度的薄层反射对褶积来合成地震记录，并针对不同噪声水平(0.1%，1%，5%，10%)，各实现500次蒙特卡洛试验，统计方法在相应噪声水平下对不同时间厚度的薄层反射可分辨的百分比(图 12b).可以看到，噪声水平和资料主频都是控制自适应滤波方法对薄层视觉分辨能力的主要因素.在高主频(宽带宽)，低噪声水平(高信噪比)下方法表现出更高的薄层分辨能力.

假设分辨比率高于50%被认为是方法可以分辨的厚度阈值，于是在不同的噪声水平下，针对不同主频的合成记录模型，统计方法可以分辨的厚度阈值，得到的结果如图 13所示.(图中黑线为不同主频的薄层调谐厚度)而按照原有算法的相同思路对加权后的方法对薄层的分辨能力进行统计的定量分析(图 13a，13b)，可以看到改进后的方法在薄层分辨能力上与原方法是相当的甚至略高于前者，这主要由于改进算法的处理结果更加稳定，这样在统计可分辨的薄层厚度阈值时定量表现更佳.

5 模型正演试验

在定量分析自适应滤波方法在分辨率和信噪比方面的应用效果之后，为了考察其对于不同地质构造的可行性及保真度，文中分别设计断层模型和河道模型进行数值模拟分析(图 14).试验以波动方程正演为手段，基于爆炸反射界面原理近似获得对应的偏移后的剖面(25 Hz子波)，并应用算法(WAPES)进行处理，同时以对应的高频子波(32 Hz)正演剖面进行对比分析.

对于断层模型，如图 15所示，自适应滤波方法的处理结果比原模拟剖面(25 Hz)表现出更细节的波形和相位信息，断层更加清晰，并且一些薄层被揭示出来，与高频模拟剖面(32 Hz)的构造特征基本一致.需要注意的是对于大倾角的构造位置处，方法的处理效果会因为子波的畸变效应变差.对于河道砂体模型，如图 16所示，在25 Hz的模拟剖面上是难以分辨五套砂体的，而处理后的剖面则表现出更丰富的信息来识别砂体，与对应的32 Hz子波模拟剖面相一致.相对于五套砂体的右侧，左侧的处理结果表现的分辨率效果要差很多，这也正是由于应用方法应用过程中选取统一子波进行处理，而高陡界面构造会造成子波的畸变引起的.因此，在方法的应用中需要考虑到选取统一的子波对方法的适用性，特别是对于高陡构造适应性的影响.

6 实际资料处理

最后，将自适应方法应用于四川某页岩探区的实际资料，对其高分辨率处理的应用效果进行探讨.从处理结果可以看出(图 17a)，高分辨算法(WAPES)有效提高了三维地震资料的分辨率，可以更精细地解释层位，并且揭示更多的储层信息，比如砂体、小断层等等(图 17b).从振幅和瞬时相位的水平切片上可以很明确地看到，处理后的资料表现出更丰富的振幅和相位信息(图 17c).进一步地，分别沿目标页岩储层位置提取相干和曲率属性来预测储层的裂缝发育情况(图 17d)，并与探区水平生产井区进行对照，可以发现，对于原始资料，探区生产井周围的相干和曲率属性未体现出较多的裂缝发育行为，而在高分辨处理后的数据上所提取的相干和曲率则在井区周围表现出了更多细节信息，这一点由井区的生产资料也得到了验证.

7 结论

本文首先对基于自适应滤波的高分辨率处理方法进行了理论探讨和模型分析，并根据地震资料噪声特点提出了加权方式的改进算法(WAPES).通过统计实验，对改进前后的方法在信噪比和分辨率方面的表现进行了定量分析.实验证明，改进后的加权方法大幅度提高了处理结果的信噪比，同时有效提高了资料的分辨率和带宽，处理效果稳定，所拓高频可信，丰富了地震资料的波形和相位信息.接下来应用典型模型试验证明方法对不同地质构造的适用性，进而将改进加权方法对某一页岩气探区进行处理，并提取裂缝相关属性(相干，曲率)对目标层的裂缝发育带进行刻画.处理结果表明，新方法处理后的数据表现出更多相位信息，提取的几何属性表现出储层裂缝发育的更多细节特征.