假设地震记录的道数为n,各地震道包含的采样点数为m,那么总采样点数为n×m个元素。地震数据矩阵r可表示为:

令x=vec(r)为矩阵r的拉伸向量,即x∈R^nm; w∈R^nm是地震数据中的随机噪声,K∈R^nm×nm是褶积算子,f∈R^nm是实际地震信号,满足如下关系:

f=Kx+w (2)

在K给定的前提下,从含噪声的实际地震信号f中重建初始地震信号x,即重建地震信号可视为反褶积和去噪同时处理。若K为单位算子,则此操作仅为去噪处理。依据Rudin等(1992),Vogel和Oman(1998)研究, ,即使用重建初始地震信号x可对公式(3)进行最小化处理:

式中: 为L¹范数; Φ_fid(x, f)是噪声信号的保真项。

本文目标是从f中重建x,但此过程对噪声w过于敏感,重建x存在不稳定性。为解决此问题,需融入TV正则项(Rudin et al,1992; Wang et al,2008),即TV的离散形式:

式中: 为L²范数; D_i为局部有限差分算子; D_ix∈R²为x在采样点i处垂直方向与水平方向的一阶有限差分。地震信号的重建过程为求解式(5)的最小化(Yang et al,2009a):

即地震信号的重建模型为:

式中:μ>0为平衡二者的参数。注意:x是理论上的真实地震信号数据; x为被恢复后的地震数据。本文通过交替方向乘子法(ADMM)求解式(6)。

ADMM的基本思想需追溯到Gabay和Mercier(1976)的研究,即令θ₁(?)和θ₂(?)为凸函数,A为连续的线性算子,则最小化能量函数为:

通过引入辅助变量v,式(7)可等价转换为:

式(8)的增广Lagrangian函数为(Gabay,Mercier,1976):

式中:λ和β为增广Lagrangian求解最优化的引入参数。因此,交替方向乘子法的迭代解如下:

式(10)的交替最小化迭代思想为:首先由初始化(v^k,λ^k)获取u^k+1,再使用u^k+1和λ^k获取v^k+1,从而可得到λ^k+1。最后,利用交替方向最小化的多次优化来获得恢复的地震信号。

首先,将式(6)转化成式(8)形式,可表示为:

式中:y=(y¹; y²)∈R^2nm, y¹和y²是长度为n×m维向量, 满足((y¹)_i;(y²)_i)=y_i∈R²。利用二次罚函数法把式(11)转成无约束优化函数形式:

(12)

式中: 为罚参数。通过ADMM对式(12)进行求解(初始x=x^k,λ=λ^k),可得到:

式中:L_A(x,y,λ)为式(12)的增广Lagrangian函数:

式中:λ_i为R²空间中的向量。根据Wang等(2008)和Yang等(2009a)研究,式(13)中argmin_yL_A(x^k,y,λ^k)等价于如下形式:

式(15)通过二维收缩来求解(Wang et al,2008; Yang et al,2009a):

在固定y^k+1和λ^k的基础上,通过最小二乘法求解式(13)中argmin_xL_A(x,y^k+1,λ^k):

(17)

式中:D=(D¹; D²)∈R^2nm×nm是一阶全局有限差分算子; D¹和D²是水平方向与垂直方向的一阶全局有限差分算子; D_i是以D¹的第i行为第1行、以D²的第i行为第2行组成的两行矩阵。为了确保式(17)中系数矩阵的非奇异性,需假设N(D^TD)∩N(K^T K)=0,其中,N(?)表示矩阵的零空间。根据Ng等(1999)研究,D^TD和K^TK经二维离散快速傅立叶变换(FFT)后可进行对角化操作处理。因此,求解式(17)可通过以下3步来完成:①计算等式右侧向量,并进行傅立叶正变换; ②通过分量形式对的特征值相除,进而获得x^k+1的二维离散傅立叶变换F(x^k+1); ③对F(x^k+1)进行二维傅立叶反变换得到新的迭代x^k+1。

最后,由式(18)更新λ:

λ^k+1=λ^k-γβ(y^k+1-Dx^k+1) (18)

式中:γ>0为更新λ过程中被赋予的步长,一般取(Wang et al,2008)。

L¹TV-SSD算法过程:输入为观测地震信号f,褶积模型算子K,μ>0,β>0以及λ>0; 输出为重建后的地震信号x。初始化条件为:x=f; λ=λ₀。当算法未收敛时,即执行以下操作:对于给定的(x^k,λ^k),根据式(16)计算y^k+1,通过求解式(17)得到x^k+1,由式(18)更新λ^k+1。算法中“not converged”满足,即通过x相对的变化来终止算法,ε>0是给定的公差。对于每次迭代,求解式(17)的主要过程为一次傅立叶正变换和一次傅立叶反变换,时间复杂度为O(m²logm)。

为了对L¹TV-SSD算法进行有效性评价,将其应用于合成地震记录上进行实验,并与维纳滤波和拉普拉斯滤波结果进行比较。实验中选取参数为:γ=1.618,ε=10^-3,罚参数β=10; 为了使地震信号重建获得最好的结果,不同的地震信号数据,μ值取值不同,范围一般在(0,60]。

采用信噪比(Signal-Noise Ratio,简称SNR)来评价地震信号重建的质量,SNR定义为:

式中:x为原始地震记录; E(x)为x的均值强度; x是重建后的地震信号数据; SNR单位为dB。

图1 不同方法对褶积后的合成地震记录的处理结果
Fig.1 Results of convoluted synthetic seismic records by different methods

原始合成地震记录主频为35 Hz、1 ms采样、共150道,每道采样点数为600,地震波为雷克子波(图1a)。使用高斯核函数(标准差σ=2)与原始地震记录褶积得到数据的SNR为7.69 dB。从图1b中可以看出,箭头指向处主频降低。使用L¹TV-SSD算法对图1b数据进行处理(μ=5),箭头指向处的分辨率明显提高,SNR值提高至13.88 dB(图1c)。图1d为SNR、优化目标函数值与L¹TV-SSD算法迭代数的关系曲线,在迭代48次后SNR持续增加,优化目标函数值随迭代次数的增加一直逐渐降低,在近140次时降到最低。

图1e与f分别是使用维纳滤波与拉普拉斯滤波的处理结果,SNR分别为10.42 dB与9.13 dB。维纳滤波的处理效果提升不大,拉普拉斯滤波处理结果其主频有所提高,但多了噪声。综上,L¹TV-SSD算法获得了较好的处理效果。

对原始地震记录(图2a)进行加噪处理后进行实验。使用高斯核函数和原始地震记录进行褶积,再通过Matlab中“randn”函数增加随机噪声的数据如图2b所示。高斯核函数的标准差σ=1.5,随机噪声系数设置为0.02。使用L¹TV-SSD算法对图2b数据进行处理(μ=4),从图2c可

图2 不同方法对褶积与含噪声的合成地震记录的处理结果
Fig.2 Results of the convoluted synthetic noisy seismic records by different methods

明显看出,L¹TV-SSD算法有效地衰减了随机噪声,提高了分辨率,SNR由5.79 dB增至11.98 dB。从图2d可以看出,随着迭代次数的增加,SNR不断提升,而目标函数值不断降低。对于维纳滤波结果(图2e)和拉普拉斯滤波结果(图2f),尽管二者的同相轴较图2b清晰,但增加了较多的噪声,使得信噪比下降,这表明这两种结果极易受噪声的影响。图3为加噪声数据、维纳滤波方法和本文方法处理结果的第125道频谱对比图。从图中可以看到,L¹TV-SSD算法提高了子波的主频,拓宽了有效频带,起到了反褶积的作用,高频随机噪声同时也得到了较好的压制; 而维纳滤波方法引入了较多的随机噪声,降低了子波的主频。

对原始地震记录(图4a)增加较强的随机噪声,其数据如图4b所示,随机噪声系数为0.15,用于褶积的高斯核函数的标准差σ=1.5。由于维纳滤波和拉普拉斯滤波受噪声影响太大,因此,本文不做讨论。图4c为L¹TV-SSD算法的处理结果(μ=4),从图中明显可看出,即使在信号受到较强噪声污染时,也可以获得不错的处理结果,

SNR由-8.48 dB增至7.60 dB。同样,从图4d可以看到,随着迭代次数的增加,SNR不断提升,而优化目标函数值不断降低。

图3 加噪声数据、维纳滤波方法和本文方法处理结果的第125道频谱对比
Fig.3 Comparison between the spectrums of the 125th trace respectively obtained by the noisy data method,the Wiener filtering method and the proposed method in this paper

图4 本文方法对褶积与强噪声的合成地震记录的处理结果
Fig.4 Results of the convoluted,strong-noise,synthetic seismic records by the proposed method in this paper

图5 不同方法对实际数据的处理结果
Fig.5 Results from different methods applied to the observation data

图6 不同方法对CMP道集数据的处理结果
Fig.6 Results from different methods applied to the CMP data

为了更好地说明L¹TV-SSD算法的适用性,将其应用于辽河某地区三维地震资料处理中。图5a为一初始叠加剖面,信噪比较低; 图5b为L¹TV-SSD算法处理结果; 图5c和d分别为应用GeoEast软件的径向预测滤波(RadiPredicFilt)和随机噪声衰减(RNA)两个模块的处理结果。3种方法的对比剖面分析表明:L¹TV-SSD算法对原始叠加剖面的背景噪声具有更好的压制效果,并且对模糊的同相轴有更好的处理效果,使其连续性增强,同相轴更加清晰,层间信息更加丰富,整体的分辨率明显提高。从RadiPredicFilt和RNA两个模块来看,RNA在噪声衰减以及细节的清晰度方面要优于RadiPredicFilt。

为了进一步验证L¹TV-SSD算法的有效性,选取实际的CDP道集(图6a)进行处理,从图6a可以看出,实际资料存在很强的随机噪声,对有效信号造成了强烈的干扰。图6b～d分别为使用L¹TV-SSD算法、径向预测滤波(RadiPredicFilt)和随机噪声衰减(RNA)模块对实际的CDP道集处理的结果。从处理效果来看,3种方法均有效地压制了随机噪声,L¹TV-SSD算法有其更好的压制效果,突显出来的细节信息要优于其它两种方法。