基于自相似性和低秩先验的地震数据随机噪声压制

程文婷,方文倩,付丽华

(中国地质大学(武汉)数学与物理学院,湖北武汉430074)

地震数据随机噪声压制方法较多,主要有基于稀疏变换的方法,如Fourier变换[1]、Wavelet变换[2-3]、Curvelet变换[4]、Radon变换[5]、Shearlet变换[6]、Seislet变换[7]等;基于滤波理论的方法,如中值滤波[8-9]、f-x反褶积[10]、聚束滤波[11]等;基于波动方程的方法[12]和基于字典学习的方法[13-14]等。近年来,矩阵降秩作为一种重要的数据处理方法,成为地震信号处理领域的研究热点[15-17]。

基于降秩理论的地震数据去噪原理是:不含噪声的地震数据具有低秩结构,而随机噪声的存在会增加数据矩阵的秩。因此,地震数据去噪问题可以转化为矩阵降秩问题。奇异谱分析(singular specturm analysis,SSA)[18]是经典的降秩方法,通过对地震数据频率切片构造的Hankel矩阵进行降秩达到去噪目的。由于奇异值分解(singular value decomposition,SVD)计算时间长,OROPEZA等[19]以随机SVD代替SVD,提高了计算效率。但是,此类方法需要事先确定秩的大小,而秩的选取将直接影响数据的去噪效果。核范数最小化(nuclear norm minimization,NNM)方法是另一类低秩矩阵近似方法。核范数作为秩函数的凸松弛函数,在其基础上构建的目标函数是凸问题,可使用凸优化方法求解,如:奇异值阈值法(singular value threshold,SVT)[20]、加速近端梯度法[21]、交替乘子方向法(alternating direction method of multipliers,ADMM)[22]、不动点迭代法(fixed point continuation,FPC)[23]等。然而NNM方法是将所有奇异值的和最小化,求出的解往往只是原始秩最小化问题的次优解。为了解决该问题,HU 等[24]提出截断核范数方法,对[min(m,n)-r]个较小奇异值的和进行最小化处理,其中m,n为矩阵的尺度,r为矩阵的秩。相较于核范数,截断核范数能更好地逼近秩。

近年来,非局部自相似先验在图像去噪领域取得了巨大的应用成果,如BUADES等[25]提出的非局部均值(non-local means,NLM)算法以及DABOV等[26]提出的三维块匹配(block matching and 3D collaborative,BM3D)算法等。事实上,地震数据相邻道之间同相轴的时差变化规律相近,地震记录在纵向和横向上存在较强的相似性,自相似块构造的矩阵具有更好的低秩性能。BONAR等[27]首次利用NLM算法压制地震数据随机噪声。SHANG等[28]提出一种自适应滤波参数估计的NLM算法,提高了自相似块的搜索效率。张岩等[29]提出非局部自相似与字典学习相结合的算法,通过对基于多道相似组的地震数据构建字典和稀疏编码,提高稀疏性压制随机噪声。ZHANG等[30]将二维相似块组排列成三维顺序结构,增强自相似块的相干特性,进而利用滤波对地震数据随机噪声进行压制。

本文提出基于自相似性和低秩先验的地震数据去噪方法,通过块匹配算法搜索地震数据的自相似块,然后以“自相似块组”为单元,利用低秩约束进行随机噪声压制。TNNR比经典的核范数更加接近原始的秩函数,但由于该问题是非凸问题,因此,我们将目标函数进行转换,然后利用APGL算法进行求解。

1 方法理论

1.1 基于降秩理论的地震数据去噪模型

地震数据去噪模型可表示为:

M=X+E

(1)

式中:M和X分别表示含噪地震数据和待恢复的无噪数据;E为随机噪声。基于降秩理论的地震数据去噪模型采用如下目标函数:

(2)

式中:rank(X)为X的秩,表示X非零奇异值的个数;‖·‖F为Frobenious范数;λ为正则化参数。

秩函数是非凸非光滑函数,导致公式(2)的求解是一个NP难问题。可利用核范数近似秩函数,将公式(2)转化为凸优化问题:

(3)

1.2 基于截断核范数的低秩模型

(4)

因为公式(4)是非凸的,所以先利用下文定理1对目标函数进行转换,然后再利用APGL算法求解。

1.3 目标函数转换

定理1[24]:已知矩阵X∈Rm×n,存在矩阵A∈Rr×m,B∈Rr×n,满足AAT=Ir×r,BBT=Ir×r,对任意非负整数r(r≤min(m,n)),有:

(5)

这样:

(6)

于是:

(7)

对公式(7)进行迭代求解。令X1=M,第l次迭代时,对Xl进行SVD分解得到Al和Bl;然后固定Al和Bl,更新Xl+1。目标函数((4)式)可转化为如下优化问题:

(8)

因此,可采用凸优化方法求解。APGL算法是经典的凸优化方法,其优势为收敛速度快,对噪声具有较强的鲁棒性[21]。本文采用APGL算法进行优化求解。

1.4 相似块匹配

相似度定义为欧式距离:

(9)

式中:mi和mj分别表示目标块和搜索的相似块,di,j越小,则mi和mj越相似。由于地震数据量较大,在全局搜索相似块计算较为耗时,因此,设置大小为Q×Q的搜索窗,在搜索窗内按照公式(9)所示的相似性度量寻找相似块。相似块大小q,相似块个数H及搜索窗Q是相似块匹配中较为重要的参数。实验中可以通过多次试误调节,在计算复杂度和去噪效果之间折衷选择参数。

由相似块组成的矩阵Mi是本文需要恢复的矩阵,因此,公式(8)转化为基于相似块矩阵形式的去噪模型:

(10)

图1 自相似块匹配与去噪示意

1.5 APGL优化求解

APGL算法是一种有效且稳定的凸优化方法,可以解决如下无约束的非光滑凸问题:

(11)

式中:g为闭凸函数;f为可微凸函数。

APGL算法将公式(11)写成如下F(X)在某一定点Y处的近似形式:

(12)

对任意t>0,迭代更新X,Y,t从而优化公式(11)。在第k次迭代,通过最小化G(X,Yk)更新Xk+1:

(13)

对于目标函数(公式(10)),令:

g(X)=‖Xi‖*

则:

(14)

对矩阵X∈Rm×n进行SVD分解:

X=UΣVTΣ=diag({σi}1≤i≤min(m,n))

(15)

矩阵奇异值收缩算子Dτ定义:

(16)

由文献[20]可知,对任意τ>0,Y∈Rm×n,有:

(17)

因此,公式(14)可转化为:

(18)

接着,更新ti,k+1和Yi,k+1:

(19)

(20)

本文所提的SP-TNNR算法具体实现过程见表1。

表1 基于自相似性和低秩先验的地震数据随机噪声压制算法实现过程

2 数据实验

利用仿真数据和实际地震数据验证SP-TNNR算法的性能,并与传统的Curvelet变换、f-x反褶积、TNNR算法以及基于自相似块匹配的核范数最小化方法(SP-NNM)进行对比分析,用信噪比(RSN)评价算法性能:

(21)

式中:S与S*分别表示不含噪数据和去噪后的数据。

2.1 仿真数据实验

仿真数据实验选取的地震数据共256道,每道256个时间采样点,采样间隔为2ms。在综合考虑算法时间复杂度、计算复杂度及去噪效果等因素的基础上,SP-TNNR算法实现中,设置搜索窗Q=30,相似块大小q=8,相似块个数H=150,超参数λ=0.95。f-x反褶积频带为1～100Hz,滤波长度为12。而在TNNR方法中秩取15。在仿真和实际数据实验中Curvelet变换尺度参数均为5,阈值大小取决于去噪数据的信噪比。SP-NNM方法中搜索窗、相似块大小以及相似块个数设置均与SP-TNNR方法一致。图2a为原始仿真叠前地震数据,图2b为加入了随机噪声(均值为0,标准差为50的高斯白噪)的仿真叠前地震数据(RSN=5.1)。图3和图4分别为5种方法去噪结果及对应的残差剖面。Curvelet变换、f-x反褶积、TNNR、SP-NNM以及SP-TNNR方法去噪后的信噪比分别为23.0,20.2,12.5,26.3,27.5,耗时分别为2.8,3.4,860.1,931.7,932.3s。从去噪结果和残差剖面中可以看出,Curvelet变换去除了绝大部分噪声,但同相轴有模糊现象且同相轴附近还残留部分噪声。f-x反褶积方法去噪结果中仍然可以明显看到噪声残留。TNNR去噪结果中同相轴附近残留噪声湮没了有效信息,去噪效果较差。SP-NNM和SP-TNNR方法去噪后噪声基本消除,但SP-NNM去噪后同相轴周围边缘过于光滑,在去噪的同时去掉了部分有效信号,而SP-TNNR去噪结果更干净,去噪后同相轴十分清晰,去噪效果最佳。

图2 原始仿真叠前地震数据(a)及加噪数据(b)

图3 5种方法去噪结果a Curvelet变换; b f-x反褶积; c TNNR;d SP-NNM; e SP-TNNR

图4 应用不同方法对仿真叠前地震数据去噪后的残差剖面a Curvelet变换; b f-x反褶积; c TNNR; d SP-NNM; e SP-TNNR

图5给出了不同噪声水平下5种方法去噪结果的信噪比。可以看出,随着噪声水平的增加,5种方法的信噪比逐渐降低,而本文提出的SP-TNNR方法去噪效果最优。

图5 不同噪声水平下5种方法去噪结果的信噪比

2.2 实际数据应用

实际叠后地震数据来自网站https:∥github.com/sevenysw/MathGeo2018,如图6a所示,包含256道,每道256个时间采样点。加入随机噪声(均值为0,标准差为50的高斯白噪)的结果见图6b(RSN=9.0),可以看出加噪数据的同相轴边缘信息模糊。经过多次试误调节,实验中搜索窗Q=30,相似块大小q=9,相似块个数H=150,正则化参数λ=0.99。f-x反褶积频带为1～100Hz,滤波长度为14。TNNR方法中秩取18。图7为实际叠后地震数据去噪结果。图7a为Curvelet变换去噪结果,该方法去掉了大部分噪声,但是同相轴过于光滑,有效信号损失较大;图7b为f-x反褶积方法去噪结果,部分有效信号和噪声同时被去掉;图7c为TNNR方法去噪结果,效果极差;图7d为SP-NNM方法去噪结果,去掉大部分噪声的同时也去掉了部分有效信号;图7e 为本文提出的SP-TNNR方法去噪结果,可以看出噪声被去掉,同相轴清晰可见且连续性较好,具有较高的保真度。5种方法去噪后的信噪比依次为19.8,18.9,13.1,20.7,21.9,耗时分别为2.3,3.2,883.2,906.3,912.8s。图8为图7对应的残差剖面。可以看出,本文提出的SP-TNNR方法去噪后有效信息得到了最大程度的保留,信噪比高于另外4种算法。图9 至图11为对应的频率-波数图。可以看出,SP-TNNR方法去噪后频率比较集中,其它4种方法去噪后依然存在频散现象。

图6 实际叠后地震数据(a)及加噪数据(b)

图7 实际叠后地震数据去噪结果a Curvelet变换; b f-x反褶积; c TNNR; d SP-NNM; e SP-TNNR

图8 实际叠后地震数据去噪后的残差剖面a Curvelet变换; b f-x反褶积; c TNNR; d SP-NNM; e SP-TNNR

图9 实际叠后地震数据(a)和加噪数据(b)频率-波数图

图10 利用 Curvelet变换(a)和 f-x反褶积(b)方法去噪后的地震数据频率-波数图

图11 利用 TNNR(a)、 SP-NNM(b)和SP-TNNR(c)方法去噪后的地震数据频率-波数图

3 结论

本文在TNNR基础上引入非局部自相似先验,使用相似块匹配方法构造低秩矩阵,提出了基于自相似性和低秩先验的地震数据随机噪声压制方法,即SP-TNNR算法。与经典的核范数方法相比,本文方法优点为:TNNR算法保留最大的几个奇异值仅惩罚较小奇异值,可以更加准确地逼近秩;与截断核范数方法相比,本文方法优点为:相似块匹配方法构造低秩矩阵充分利用地震数据的自相似性,通过相似块构建的矩阵具有更好的低秩性能,加入自相似先验的截断核范数方法比没有增添自相似先验的截断核范数方法具有更好的随机噪声压制能力;与经典的Curvelet变换和f-x反褶积方法相比,本文方法优点为:SP-TNNR方法去噪后的数据信噪比更高,去噪后的地震数据纹理细节更加完整,有效信号充分保留。

仿真和实际地震数据实验结果均表明本文提出的SP-TNNR方法相较于其它方法去噪效果更好,具有更高的信噪比。但本文方法也存在着不足,相似块搜索耗时较长,导致本文方法计算效率较低。快速的搜索算法以及相似块匹配中参数的科学选择方法将是未来的研究方向;另外,TNNR算法能更加逼近秩,但求解过程中需要两层SVD分解,计算复杂度较高,如何提高计算效率也是下一步的工作。

致谢:本文在完成过程中采用了哈尔滨工业大学数学学院马坚伟教授团队的MathGeo2018工具包,在此表示感谢!