信号的稀疏表示至关重要,离散余弦变换、小波变换、曲波变换等都能稀疏表征信号,因此常常被用作压缩感知的稀疏变换基。曲波变换具有各向异性和多尺度、多方向的诸多优点,较其它2种方法稀疏效果更好(Herrmann,Hennenfent,2008),是地震数据稀疏表达最有效的方法之一,用离散曲波变换来稀疏表征地震信号的方法原理如下:

假定输入信号为f [t₁,t₂],其中t₁≥0,t₂<n,t₁<t₂,基于wrapping算法核心思想:围绕原点包裹,对任意区域映射到原点的仿射区域(杨连刚等,2019),其具体变换流程如下:

①将f [t₁,t₂]通过二维的傅立叶变换和FFT域采样得到表达式为:

②每个尺度j和角度l与窗函数作用后得到:

③对于每一个大小为L₁×L₂的窗数据,围绕其原点,可得:

④对每个做二维傅立叶逆变换得到曲波域的系数表示为:

式中:表示离散波形。

⑤取软硬阈值如下:

式中:σ是阈值。

压缩感知理论核心思想是如果信号是可压缩的或在某个变换域是稀疏的,用一个与变换基不相关的观测矩阵将高维信号投影到低维空间上,然后通过重构算法就可以从低维信号重构出原始高维信号(Candes et al,2005)。

假设气枪信号为y(n),其重构问题可以通过以下的模型进行表示:

y(n)=Rvec[x(n)]　(6)

式中:y(n)表示通过测量矩阵得到的不完整的地震数据; vec表示二维的向量化算子; x(n)表示二维完整地震数据; R∈R_PXMN(P≤MN),是根据缺失数据定义的单位测量随机矩阵,假设O为缺失数据位置,R=0; I为未缺失位置数据的位置,R=1(Kong,Yu,2012)。在定义测量矩阵时,R必须是随机测量矩阵,根据压缩感知理论,如果x-在某一个变换域Ψ中是稀疏的,其变换系数为:

式中:α是Ψ的等价或者稀疏表示。

根据缺失的位置设计一个平稳的、变换基Ψ不相关的M×N维的测量矩阵[WTHX],这里选用一个服从Z(0,1/Z)、正态分布、大小为M×N高斯随机矩阵来满足压缩感知理论对于随机投影矩阵的要求。对α进行观测得到观测数据为:

为欠定方程组,有多个方程解,但α具有稀疏性,可以通过不断促进α的稀疏性进行反演求解。这就为实现缺失数据的精确重构成为可能。其中的稀疏促进策略(Herrmann,Hennenfent,2008)如下:

式中:是α的估计值,λ是阈值参数。

图1 L1范数(a)和L2范数(b)求解示意图(p表示范数)
Fig.1 Schematic diagram of solving L1 norm(a)and L2 norm(b)

模型空间限制在x的一个L-ball 中。在(x₁,x₂)平面上可以画出目标函数直线,而约束条件则成为平面上大小不断变化的一个norm ball 。直线与 norm ball 首次相交的地方就是最优解。与L₂范数相比,L₁范数更有可能得到值为0的解(图1),使得x更具稀疏性(姚华建,2013)。

本文利用L₁范数求解得到,通过不断促进α的稀疏性来得到的最佳估计值。是α的L₁范数,[HJ0]利用“0-范数”意义下的优化问题求解α的[HJ2.2mm]精确或者近似逼近。通过式(9)估计出最佳系数值以后,可以通过α=求出缺失地震数据y(n)的最佳重构数据。

气枪信号可以表征为y(n)=r+n,其中,r代表有效信号,n表示为噪声信号。基于曲波域的阈值滤波方法(CSCT),域中小系数容易被过度去除,缺失信号关键细节信息,并且不连续的区域常会出现震荡的伪影现象(唐刚,杨慧朱,2010),影响数据的降噪资料。因此,需要首先利用压缩感知重建策略对因阈值滤波导致缺失的系数进行重建,解决信号稀疏特征的连续性问题,再通过全变差最小化技术对重建后的系数进行调整(Newark,1998; 沈维燕,2006),达到联合降噪的目的。

利用全变差最小技术对重建的曲波系数进行调整主要有两种方式:一种是调整大于设定阈值的系数:

c^D_n+1=c^D_n-tⁿp_≥(g_TV(c^D))　(10)

另一种是调整小于阈值设定的系数:

c^D_n+1=c^D_n-tⁿp_&lt;(g_TV(c^D))　(11)

式中:p_&lt;表示在小于设定阈值的曲波系数空间上的投影; p_≥表示在大于设定阈值的曲波系数空间上的投影; 设置最大迭代次数N和初始值,n=1。g_TV(c^D)通过以下过程求解。

假设Ω是R^N中的有限子集,曲波域中系数c^D的全变差定义为:

其中

对上式进行分部积分得到系数函数的全变差为:

式中:Ω为图像c^D的支撑区间; x∈Ω为图像的坐标向量。

基于全变差的降噪方法可通过最小化E(c^D)来实现:

式中:λ是拉格朗日常数,控制惩罚项的权重; 目标函数E(c^D)是c^D的凸函数,其存在的一阶充分必要条件是ΔE(c^D)=0,因此,其对应的拉格朗日方程是

假设其满足Neumann边界条件,可通过投影梯度算法进行求解:

式中:初始值可设为c^{D_[HT4.5]1[HT]}=c^D₀; t≥0表示迭代步长; g_TV(c^D)表示全变差函数的次梯度。

式中:β>0为一个很小的正值,以解决全变差函数在某些点的不可微性,取步长

式中:n=n+1,判断n是否接近N,如果是迭代结束,否则继续求解次梯度。

基于曲波城的压缩感知降噪方法具体的流程如图2所示。

图2 基于曲波域的压缩感知降噪方法流程图
Fig.2 The compress sensing denoising method based on curve wave domain

基于曲波域全变差的压缩感知降噪方法对低信噪比的气枪信号实现了高效稳定的降噪效果。为了验证该方法应用于气枪信号降噪的处理效果,其中,已知的理论一维信号长度为1 400个采样点,信号带宽2～8 Hz,采样率100 Hz,且具有两个正弦周期(向涯等,2019)。在每道信号上随机加入6 dB高斯白噪声,再从100道中随机选择其中的9道按照线性递增的方式再次加入不同比例的随机噪声,噪声最高信号信噪比达到-15 dB,如图3a所示,从图中可以看出,加入噪声的9道信号相对原始信号已经基本看不出有效信号。

按照本文方法的处理流程,对加噪声信号进行降噪处理,在参数的选择过程中,其对曲波域阈值的选择综合对比了3种阈值(硬阈值、软阈值、软硬阈值结合)的降噪效果,其中软硬阈值结合的降噪效果较好,折衷系数选择0.001,能够最大限度保留信号的细节信息。在对曲波域系数进行压缩重建过程中,信号的缺失率不能低于30%,否则信号不能得到有效恢复。迭代次数选择也尤为重要,次数较大会引入随机噪声,较小对信号的细节信号恢复效果较差,综合尝试后认为选择迭代6次较为合适。处理结果如图3b所示,从图中可以看到,即使信噪比很低的信号,随机噪声都能得到很好地压制,与模拟信号的波形和振幅一致性仍较强,相位恢复精度较高,获得较好的处理效果。

图3 模拟信号加载不同比例噪声(a)以及本文方法降噪结果(b)
Fig.3 The original signal loaded with different proportions of noise(a),and the noise reduction by the method in this paper(b)

为了进一步验证该方法的可行性和优势,本文引入带通滤波(BDF)、小波变换滤波(DWTF)、S变换模板滤波(STTF),分别对加噪声信号进行处理,过程如下:用带通滤波对信号进行2～8 Hz带通滤波。小波变换滤波利用小波函数对信号进行小波3层分解后,其中阈值统一取1,进行阈值滤波。S变换模板滤波对信号进行叠加后,变换到时频域中,通过设定时频域区间的阈值窗口得到滤波模板,利用滤波模板再对信号进行滤波。4种方法的处理结果如图4所示,选择其中一道进行处理效果评价(图5)。

如图5所示,通过对4种方法处理结果的信噪比SNR和均方根误差MSE可以看出,随着噪声比例越高,传统的几种方法处理后的信号SNR越低,均方根误差MSE越高。而本文的方法对不同比例的噪声信号降噪效果明显高于其它方法,即使信噪比很低也能达到很好的效果,并且稳定性更强。从单一道的处理效果(表1)来看,BDF方法可以将高频成分去除,但是与有效信号频率相近的噪声依然存在。用DWTF方法处理后虽然有效信号得到加强,但是高频成分随机噪声依然存在,STTF方法虽然克服了前两者的缺点,但是部分相位没有得到有效恢复,造成细节缺失,而本文的方法不仅很好地压制了随机噪声,而且恢复了对相位细节。

图4 不同滤波方法与本文方法处理效果对比
Fig.4 Comparison of processing effects between various filtering methods and the method in this paper

表1 不同滤波降噪方法的处理结果对比评价
Tab.1 Comparison and evaluation of processing results of different filtering and noise reduction methods

图5 不同滤波降噪方法处理结果的信噪比(a)和均方根误差对比(b)
Fig.5 Comparison of the results of SNR(a)and MSE(b)by different filtering and noise reduction methods

为了进一步验证本文方法在波速测量中的高效稳定性,由已知一维气枪信号利用移动窗口压缩-拉伸法(Meier et al,2010; 刘志坤,2010)生成波速变化率不超过±5%的100道二维模拟信号(图6a)。在每道信号上随机加入6 dB高斯白噪声,再随机从100道中选择其中的9道按照线性递增的方式再次加入不同比例的随机噪声(同2.1节),如图6b所示。按照上述的处理流程分别对加噪声信号进行降噪处理,对含噪信号进行2～8 Hz带通滤波。利用小波函数对信号进行小波3层分解后,阈值统一取2,进行小波阈值滤波。对信号所有道进行叠加后,变换到时频域中,通过设定时频域区间的阈值窗口,得到滤波模板,从而达到利用叠加信号进行模板滤波的目的。

图6 使用不同滤波降噪方法对二维波速变化率的模拟信号处理后的降噪结果对比
Fig.6 Comparison of noise reduction results of analog signals with two-dimensional wave velocity change rate by different filtering and noise reduction methods

由图6c可以看出,本文方法不仅可以很好地压制随机噪声,也去除了随机的线性干扰; DWTF虽然恢复了有效信号的大部分能量,但是对高频成分的噪声压制力较弱(图6d); STTF方法对随机噪声的压制力较强,但是由于模板阈值的选取存在一定的误差,导致对信号细节的恢复较差,并且对加入的较强随机线性干扰的压制能力较弱,一些线性干扰依然存在(图6e); BDF方法只能对固定频段的噪声进行压制,高频噪声和线性干扰依然存在(图6f)。

为进一步验证本文方法的稳定性和高效性,对上述几种降噪方法处理后的信号进行波速测量精度评价,如图7所示。

从图7可以看到,利用本文方法降噪后测量的信号波速变化率与模拟信号(蓝线)的量级一致,波速变化曲线基本一致,且拟合度较高。用DWTF方法处理后测量的信号波速变化率虽然与模拟信号的趋势一致,但是依然存在一定的误差,并且在较强的线性干扰处存在突变,用BDF方法处理后测量的信号波速变化率由于对随机噪声和线性干扰压制较弱,导致波速测量的精度严重不足,并且在较强的线性干扰处存在突变量。

图7 不同滤波降噪方法处理后波速测量结果
Fig.7 Wave speed measurement results after processing by different filtering and noise reduction methods