石战战，庞溯，王元君，池跃龙，周强

(1.乐山师范学院 人工智能学院，四川 乐山 614000；2.成都理工大学 工程技术学院，四川 乐山 614000；3.成都理工大学 地球物理学院，四川 成都 610059；4.西华师范大学 国土资源学院，四川 南充 637002)

0 引言

地震信号不可避免受噪声干扰[1-3],有效衰减地震噪声能够为后续处理流程提供高保真叠前道集，能够提高速度分析、偏移和反演等处理环节的准确度。地震噪声包括随机噪声和相干噪声2类，不同类型地震噪声具有不同的传播规律和分布特征，需要采用不同的方法或算法进行压制[4]。其中,随机噪声表现为随时间变化的高频随机震荡，会降低地震信号的分辨率和同相轴的连续性[5]。

随机噪声压制是地震信号处理的常规环节之一，虽然业内已经提出了众多的去噪算法，但同时实现最大程度压制随机噪声和最小程度损害有效信号依然是地球物理学界的一个开放问题[4]。已有的随机噪声压制算法可以归纳为5类：①预测滤波算法，利用了信号空间方向具有可预测性，主要包括f-x域反褶积[6]、f-xEMD预测滤波[7]、t-x域预测滤波[8]等方法。这类算法的优点是同时处理多道信号，能够有效利用信号的道间相干性，有效压制随机噪声，缺点是容易损害有效信号，尤其是陡倾同相轴。②基于信号分解的去噪算法，将含噪声信号分解为不同分量，舍弃噪声能量占主导的分量实现噪声压制；常用的分解算法主要有经验模态分解(empirical mode decomposition, EMD)[7,9-10]和变分模态分解(variational mode decomposition,VMD)[11]2类。③基于稀疏变换的去噪算法，基于地震信号可以在变换域中被压缩表示，有效信号表现为高幅度系数，而随机噪声为低幅度系数，采用阈值函数就能实现噪声衰减；常用的变换算法包括Fourier和短时Fourier变换[12]、S和广义S变换[13]、wavelet变换[14]、curvelet变换[15]和稀疏表示[16]等。④基于低秩方法的去噪算法，这类方法对地震信号作时频分析或Hankel变换后，假设纯净完整的地震数据具有低秩特征，而混叠噪声和缺失地震道会使数据的秩增加;主要包括矩阵低秩方法[4,12,17]和高维张量低秩去噪[18-19]2类,其缺点是矩阵低秩方法本质是单道去噪不能有效利用信号的道间相干性，而高维去噪算法计算复杂度过高。⑤基于深度学习的去噪方法，利用卷积神经网络(convolutional neural network, CNN)[20-23]、自编程机(Autoencoders)[24]或生成对抗网络(generative adversarial networks, GAN)[25]等网络结构实现随机噪声压制,但深度神经网络依赖大量带标签的训练样本，限制了深度学习去噪算法的广泛应用。

结合多类方法构造合成去噪算法是近年来随机噪声压制领域的热点研究方向，尤其是稀疏变换与f-x域算法和低秩、降秩算法的结合，其目的是结合2类算法的优势，克服各自算法缺陷或假设条件。合成去噪算法主要包括：①结合f-x域算法与信号稀疏变换，如混合时频分析技术[13]能够同时压制随机噪声和陡倾相干噪声；②结合f-x域算法与信号分解算法,如f-xEMD[9]、f-xVMD[11]和f-xEMD预测滤波[7]克服了f-x域算法的平稳信号假设条件；③结合f-x域算法与低秩算法，如Hankel低秩近似[26]能够有效利用地震信号的f-x域特征，提高算法去噪能力；④结合信号分解和低秩算法，如经验低秩近似算法[27]通过EMD将地震信号分解为若干低秩分量，避免了低秩近似算法的小分析窗要求；⑤结合稀疏分解和低秩算法，如地震信号的时频谱具有低秩特征，时频域低秩矩阵近似[4,12,17]能够结合时频分析和低秩矩阵近似2种算法优势，提高去噪能力。

低秩矩阵算法主要包括低秩稀疏分解(low-rank and sparse decomposition)和低秩近似(low-rank matrix approximation)。其中，低秩矩阵近似一般采用核范数正则化项，利用奇异值阈值(singular value thresholding, SVT)实现快速求解[28]。但SVT不能准确估算矩阵的非零奇异值，而采用非凸正则项能够提高奇异值估算的准确性[29-31]。常用的非凸惩罚函数包括加权核范数[32]、Transformed Schatten-1[33]和p范数[31]3种。引入非凸正则化的缺点是低秩矩阵近似问题的目标函数变为非凸函数，难以优化，且受局部最小值干扰。Ankit和Selesnick[29]提出合理构造非凸正则项，同时保证目标函数为严格的凸函数。

共偏移距道集具有平缓(甚至近似水平)的同相轴结构，满足f-x域去噪线性同相轴假设前提。结合f-x域去噪、时频分析和非凸正则化低秩矩阵近似这3种算法的优势，提出基于f-x域时频低秩矩阵近似的共偏移距道集去噪方法。该方法能够同时克服f-x域去噪的平稳信号假设和时频域稀疏低秩去噪算法难以利用信号的道间相干性这2种缺陷，同时利用多元回归理论对噪声标准差进行估计，估算出最佳的奇异值阈值。通过数值模拟和实际地震资料试算验证，所提方法能够有效压制地震随机噪声，同时保持有效信号不受损害。

1 方法原理

1.1 共偏移距道集f-x域去噪

假设地震数据中混叠随机噪声，将有效信号和随机噪声分别记为s(t,x)和n(x,t)，则叠前地震道集可表示为

y(t,x)=s(t,x)+n(t,x),

(1)

式中:t和x分别表示时间坐标和偏移距。进一步假设地震信号由单个线性同相轴构成：

y(t,x)=w(t-x/V)+n(t,x) ,

(2)

式中：w(t)为地震子波，V为速度。对式(2)等号两端作Fourier变换，并利用Fourier变换的线性性质可以得到：

Y(f,x)=W(f)eiωx/V+N(f,x) ,

(3)

由式(1)～(3)可知，t-x域线性同相轴映射为f-x域中的复谐波，且t-x域中地震信号为有效信号s和随机噪声n的叠加，这种加性关系在f-x域仍然成立。因此，对于t-x域中单一线性同相轴模型，f-x域中有效信号表现为规则谐波震荡，而随机噪声则被映射为高波数随机震荡。通过将信号变换到f-x域，能够增加信号和噪声的分布差异。类似地，当地震信号由p个线性同相轴组成时，由Fourier变换的线性性质可知，f-x域中地震信号表现为p个谐波的叠加，因此，对每一单频分量做滤波处理就能够分离地震噪声。

传统的f-x域去噪处理共炮点道集或共中心点(common mid-point, CMP)道集，而这2种道集中反射波为双曲线同相轴，直达波、折射波和面波为倾斜同相轴，倾角各不相同，且面波存在频散特征。因此，共炮点道集和CMP道集均不满足f-x域去噪假设前提，直接采用f-x域去噪会损害倾斜(尤其是陡倾)同相轴。而共偏移距道集具有平缓甚至接近水平的同相轴结构，基本满足f-x域去噪的线性同相轴假设前提。因此，在共偏移距道集中执行f-x域去噪，其去噪能力和有效信号保持能力必然会提高。

1.2 f-x域时频低秩矩阵近似去噪方法

上述推导过程建立在固定子波和线性同相轴基础上，假设地震信号为稳定信号。但实际介质为非均匀各向异性粘弹性介质，地震子波随能量传播而变化，因而，地震信号不能被视为平稳信号，要求引入非平稳信号分析方法处理每一单频分量。

时频分析能够从时间和频率2个维度对非平稳信号的变化特征进行度量。常用的时频分析方法主要有短时Fourier变换、广义S变换、wavelet变换和Wigner-Ville分布等。其中，Wigner-Ville分布为chen类分布，虽具有高时频分辨率的优点，但受交叉项干扰，在地震信号分析处理领域应用较少。短时Fourier变换、广义S变换和wavelet变换为线性时频分析方法，受测不准原理制约，不能同时实现最佳的时间分辨率和频率分辨率，但这类方法不受交叉项干扰，并且具有较高的计算效率，在地球物理学领域广泛应用。其中尤以广义S变换同时兼具较好的时频分辨率，且可以借助快速Fourier变换高效实现，本研究中将其用于分解f-x域中的每一单频分量。

由于f-x域中地震信号为谐波的叠加，因此，记f-x域中某一单频分量为Y(f0,x)，当Y(f0,x)不受噪声干扰时，其广义S变换系数矩阵具有稀疏低秩特征，而当Y(f0,x)受随机噪声干扰时，随机噪声会使广义S变换系数矩阵的秩增加。因此，求解广义S变换系数矩阵的低秩近似矩阵就能实现随机噪声分离。

f-x域时频低秩矩阵近似去噪方法的计算流程为：①由共炮点道集或CMP道集抽取共偏移距道集；②逐道作Fourier变换，将共偏移距道集变换到f-x域；③对f-x域中每一单频分量Y(f0,x)作广义S变换；④对广义S变换系数矩阵作低秩矩阵近似，分离随机噪声；⑤作逆广义S变换重构纯净单频分量Y*(f0,x)；⑥逐道作逆Fourier变换，将共偏移距道集变换到t-x域。步骤③对f-x域中每一单频分量的广义S变换系数矩阵作低秩矩阵近似计算，能够利用前述SVT优化求解，无需将地震数据变换到高维空间，可提高去噪算法的计算效率。

1.3 非凸正则化低秩矩阵近似

低秩矩阵近似的目的是从含噪声矩阵Y中估计纯净的低秩矩阵X:

Y=X+N, (X,Y,N∈Rm×n)

(4)

式中:N为噪声矩阵，R为实数集，m和n分别为矩阵的行、列数。低秩矩阵近似是一个欠定的逆问题，需要引入正则化项φ(·)对解空间进行约束：

(5)

当φ(x)=|x|时，式(5)为标准的核范数最小化问题，可以通过SVT优化求解。如前所述，标准核范数最小化问题不能准确估算非零奇异值。一种有效策略是适当选择非凸正则化项φ(·)，同时保证目标函数Ψ(X)为严格凸函数。所选非凸正则化项φ(·)应满足文献[29]中Assumption 1，其中一种常用的形式为

(6)

φ(·)的临近算子(proximal operator)定义为

(7)

当按照式(6)、式(7)取值时，最优化问题(式(5))存在唯一的全局最优解：

(8)

式中：Y=UΣVT为Y的奇异值分解(singular value decomposition, SVD);a取0≤a<1/λ。

1.4 噪声估计和最佳奇异值阈值

噪声标准差σ可以通过多元回归理论(multiple regression theory)进行估计。定义yi=[Y]:,i([Y]:,i为Y的第i列)和Y∂i=[y1,…,yi-1,yi,yi+1,…,yn]。假定yi可以由Y其余各列的线性组合表示：

yi=Y∂iβi+ξi,

(9)

式中:ξi为误差项。对每一个i∈{1,…,m}，回归向量βi的最小二乘估计可以表示为

(10)

估算噪声为

(11)

2 数值模拟

2.1 核范数正则化与非凸正则化对比

通过数值模拟对比核范数正则化与非凸正则化。实验首先生成2个零均值单位方差的Gaussian独立同分布(i.i.d.Gaussian distribution)随机矩阵L100×10和R10×100，然后生成低秩矩阵Y0=LRT。向Y0中混入σ=6的随机噪声，对比传统低秩矩阵近似(采用核范数正则化)和非凸正则化低秩矩阵近似(采用非凸正则化函数)计算奇异值(图1)。图1中绿色点划线和黑色虚线分别表示纯净低秩矩阵和噪声矩阵的奇异值，可见当矩阵混叠随机噪声后其秩必然增加;图中蓝色虚线和红色实线分别为核范数正则化和非凸正则化低秩矩阵近似2种算法估算的奇异值，对比可见非凸正则化算法估算的奇异值与真实奇异值更为接近，估计出矩阵的秩与真实值相同。

图1 核范数正则化与非凸正则化对比分析

2.2 地质模型与地震正演

为了验证所提基于f-x域时频非凸正则化低秩矩阵近似的共偏移距道集去噪方法的有效性，建立如图2a所示的5层地质模型,地层速度分别为2 000、2 300、2 600、2 900和3 200 m/s，地层埋深分别为200、500、800和1 100 m。实验采用中间放炮、双边接收观测系统模拟地震采集，道间距和炮间距均为50 m，每炮48道接收，共激发48炮。图2b、c分别为第1炮正演共炮点道集和加噪炮集。为了增加实验难度，逐道加入不同幅度的随机噪声，含噪声地震道信噪比在1～10 dB范围随机变化。共炮点道集中反射波均为双曲线型同相轴，不同深度反射波倾角各不相同，因此，共炮点道集不能视为由线性同相轴叠加而成，不满足f-x域去噪假设前提。图2d为偏移距为50 m的共偏移距道集，道集中各地层反射波均为水平同相轴，满足f-x域去噪假设前提。因此，在共偏移距道集中执行f-x域去噪，其结果的精度必然提高。

图2 地质模型及其正演数据

2.3 共偏移距道集去噪算法对比

图3所示为地震道集抽取共偏移距道集，在共偏移距道集中试算f-x域时频非凸正则化低秩矩阵近似、时频域低秩矩阵近似和f-x域反褶积3种算法。f-x域时频非凸正则化低秩矩阵近似采用1.2节所述算法流程，去噪结果和滤波残差分别如图3a、b所示。与图2d所示的含噪共偏移距道集对比可以发现，处理后随机噪声得到压制，剖面信噪比提高，无明显残留随机噪声(图3a)，残差剖面符合随机噪声分布规律，无连续同相轴泄露(图3b)，说明f-x域时频非凸正则化低秩矩阵近似算法兼具良好的去噪能力和有效信号保持能力。图3c、d分别为时频域低秩矩阵近似算法去噪结果和滤波残差，该算法逐道处理各道地震数据，对单一地震道作广义S换后，再对时频谱做低秩矩阵近似计算。可以看出，处理后部分地震道残留随机噪声(如图3c中箭头所示)，残差剖面存在连续同相轴(如图3d中箭头所示)，说明单道算法不能有效利用信号的道间相干性，去噪能力和有效信号保持能力有限。图3e、f分别为f-x域反褶积去噪结果和滤波残差，共偏移距道集满足f-x域去噪的线性同相轴假设条件，该域中计算f-x域反褶积能够同时实现压制随机噪声和保持有效信号不被损害，残差剖面上无残留连续同相轴(图3f)，但该算法去噪不彻底，处理结果中仍残留小幅度随机噪声(图3e)。通过对比可知，所提出的f-x域时频非凸正则化低秩矩阵近似具有明显优势。

图3 共偏移距道集去噪算法对比

抽取图3所示去噪结果的第4道地震数据进行对比分析(图4)。可以看出：f-x域时频非凸正则化低秩矩阵近似具有较好的保幅去噪能力，处理结果无明显残留噪声(图4c)，残差剖面无有效信号泄露(图4d)；时频域低秩矩阵近似去噪算法会损害有效信号，图4e中箭头所示反射波幅度低于真实幅度(图4a)，残差剖面泄露有效信号(图4f中箭头所示)；f-x域反褶积压制随机噪声能力有限，处理后仍残留较大幅度随机噪声(图4g)。为了定量对比3种算法的有效性，计算了3种算法去噪结果的相对误差，分别为11.697%、37.418%和21.913%，进一步说明了本文算法具有明显优势。

图4 单道对比3种共偏移距道集去噪算法

2.4 共偏移距道集和共炮点道集处理结果对比

将图3a所示共偏移距道集f-x域时频非凸正则化低秩矩阵近似算法处理结果反变换回共炮点道集。在共炮点道集对比共偏移距道集算法处理结果和滤波残差如图5a、b所示，同样可以说明共偏移距道集近似满足f-x域去噪算法的线性同相轴假设，处理后随机噪声得到压制，且有效信号保持良好。为了进一步说明共偏移距道集算法的优势，在共炮点道集中对f-x域时频非凸正则化低秩矩阵近似算法进行试算。去噪结果和滤波残差分别如图5c、d所示，可以发现，时频稀疏低秩近似算法虽然具有较好的去噪能力，但因共炮点道集不满足f-x域去噪算法假设前提，两者结合不能发挥各自算法优势，处理后剖面仍残留随机噪声(图5c)，有效信号保持能力不足，残差剖面泄露连续同相轴。对比结果说明，共偏移距道集近似满足f-x域去噪算法的线性同相轴假设，在该域中进行去噪处理能够提高算法的处理效果。

图5 共偏移距道集和共炮点道集处理结果对比

3 应用实例

利用某区实际地震数据，验证所提f-x域时频非凸正则化低秩矩阵近似算法的有效性。该数据震源为30.06来复枪，井深25 cm；测量仪器采用L40A 100 Hz检波器和BISON 24096型96道工程地震仪。图6a为第1001炮共炮点道集，可以看出剖面受随机噪声干扰严重(图中箭头所示)，同时，该数据还受面波和低频干扰污染，共炮点道集中，直达波和面波具有倾斜同相轴特征，反射波为双曲线同相轴，不满足f-x域去噪的线性同相轴假设前提。图6b为偏移距为15 m的共偏移距道集，剖面表现为平缓的同相轴结构特征，各波型均为近似水平的同相轴。因此，共偏移距道集基本满足f-x域去噪的线性同相轴假设前提，在该域中进行去噪处理，能够发挥去噪算法的优势。

图6 实际地震数据

抽取共偏移距道集试算3种算法(图7)。图7a、b分别为f-x域时频非凸正则化低秩矩阵近似算法去噪结果和滤波残差，与原始共偏移距道集(图6b)对比可以发现，该算法能够有效压制随机噪声和低频干扰，去噪后剖面信噪比提高，无残留随机噪声(图7a)；且该算法在去噪同时未损害有效信号，残差剖面无连续同相轴泄露(图7b)，说明该算法能够结合f-x域去噪与时频域低秩矩阵近似2种算法的优势，提高去噪和有效信号保持能力。图7c、d为时频域低秩矩阵近似算法去噪结果和滤波残差，结果显示随机噪声得到一定压制，但剖面仍残留小幅度随机噪声(如图7c中箭头所示)，残差剖面符合随机噪声分布规律(图7d)，其原因是传统时频域低秩矩阵近似算法处理单道地震信号，没有利用信号道间相干性。图7e、f分别为f-x域反褶积去噪结果和滤波残差，该算法能够有效压制随机噪声，但有效信号保持能力不足，处理后同相轴结构发生改变(如图7e中箭头所示)，其原因是f-x域反褶积要求待处理信号为平稳信号，而实际地震信号为非平稳信号。通过实际数据试算对比，进一步验证了所提算法能够有效压制随机噪声，同时保持有效信号不被破坏。

图7 实际数据共偏移距道集去噪算法对比

为了突出共偏移距道集地震去噪算法的优势，分别在共偏移距道集和共炮点道集试算本文算法。为了便于对比，将共偏移距道集算法处理结果和滤波残差变换到共炮点道集，图8a、b显示处理后随机噪声得到压制，剖面无大幅度随机噪声残留，残差剖面无明显连续同相轴泄露。图8c、d为共炮点道集处理结果和滤波残差，剖面残留随机噪声幅度大于共偏移距道集算法，残差剖面残留明显的连续同相轴。以上结果说明共偏移距道集同相轴结构平缓，甚至接近于水平同相轴，较共炮点道集和CMP道集相比更为符合f-x域去噪假设前提，在共偏移距道集进行去噪处理能够较好地发挥算法性能。

图8 实际数据共偏移距道集和共炮点道集处理结果对比

4 结论

1)f-x域去噪算法假设地震剖面由线性同相轴叠加而成，共炮点道集和CMP道集中直达波、折射波和面波为倾斜同相轴，反射波为双曲线同相轴，且不同地层各波型同相轴倾角不同，甚至互相干涉叠加。地震道集不能视为由线性同相轴叠加而成。当处理共炮点道集和CMP道集时，f-x域去噪算法会损害倾斜(尤其是陡倾)同相轴。而共偏移距道集具有平缓甚至接近水平的同相轴结构，基本满足线性同相轴假设条件，在共偏移距道集道集中执行f-x域去噪能够提高算法性能。

2)f-x反褶积等f-x域去噪算法假设地震信号为平稳信号，而实际地震信号为非平稳信号，不满足算法假设条件，需要将非平稳信号去噪算法引入f-x域。时频分析是一种有效的非平稳信号分析方法，且时频域数据具有稀疏低秩特征。将时频稀疏低秩近似去噪引入f-x域，能够克服传统f-x域去噪算法的平稳信号假设条件。多道地震信号同时处理，能够利用信号道间相干性提高去噪能力。