0 引言

目前地震资料随机噪声常规压制方法主要包括空间域压制和变换域压制。前者有中值滤波、各向异性扩散滤波等；后者有小波阈值法、曲波变换阈值法等。其本质都是利用在特定域中有效信号与噪声具有可区分的特征进行信噪分离。

时频分析类方法主要根据时频域中有效信号与噪声的能量分布，进行噪声的分析和定位及剔除或压制^[1-3]，主要有基于短时傅里叶变换、连续小波变换和S变换等方法^[4]。Morlet^[5]提出的改进小波变换可在时频域中寻找目标信号的局部特征。但基于小波变换的一系列延伸方法对具有复杂同相轴的地震资料，效果并不理想。针对该问题，Candès等^[6]提出了具有多方向、多尺度的曲波变换。包乾宗等^[7]将曲波变换与奇异值分解方法相结合，对地震资料中的不同成分进行分离；Neelamani等^[8]将曲波变换用于三维地震数据的随机噪声的压制；孟阁阁等^[9]基于曲波变换，设计了一种自适应阈值去噪方法，效果比较明显；张恒磊等^[10-11]采用非线性阈值法基于曲波变换对叠前资料进行了去噪处理，随后提出了适用于曲波域和空间域的倾角扫描去噪技术，并在复杂区块资料中取得了较好的效果；张之涵等^[12]提出一种基于三维曲波变换的噪声压制方法，通过改进非线性阈值，更好地保护了有效信号；姚振岸等^[13]将曲波变换与各向异性扩散滤波相结合，改善了常规曲波变换降噪方法的过度平滑与环绕效应，实现了数据边界和有效信号振幅的保护。

近几年，诸如奇异值分解、稀疏表征等新方法也逐步应用于地震资料的噪声压制。它们通过构建适用于目标资料的字典，替代小波变换、曲波变换中的确定性基函数，其本质是将去噪的过程转化为优化问题。基于该思路，Sanchis等^[14]通过独立分量分析算法构建基函数，但该方法需要有效成分的统计估计；李海山等^[15]采用形态分量分析方法进行字典设计，通过求解表征目标函数达到去噪的效果，但测试结果显示仍存在一定的伪影现象；Tang等^[16]通过引入学习型过完备字典实现最优稀疏表征，并采用字典训练与总变差量的最小化共同抑制伪吉布斯效应，避免了伪影的干扰，但该方法需要已知噪声的统计信息，在实际应用中需要前期统计与实验；徐彦凯等^[17]提出了基于自适应奇异值分解的局部正交压制方法，通过将残余信号从随机噪声中的二次分离，提高了资料的信噪比。

以上时频分析方法受不确定准则限制，只能在时间分辨率与频率分辨率之间采取折中处理^[18]。针对这一缺陷，有学者提出了用于处理时序非稳态信号的经验模态分解(EMD)方法^[19-20]。基于该方法的时频分析可以在时间、频率尺度上均达到较高的分辨率。EMD方法在地震资料噪声压制和属性分析中均获得了较好效果^[21-23]。而对于复杂的时序非平稳信号，仅在时间域进行分解往往效果不佳，或出现假频等现象，而采用EMD方法，将单道地震记录或剖面在合适的域内分解为一系列本征模态分量(IMF)，可以表征不同频段内的局部能量异常。传统EMD方法的适用条件可参见文献[21-24]。Bekara等^[23]率先将EMD方法用于频率—空间域中的各频率切片数据，将分解得到的第一模态作为噪声剔除，重构之后得到噪声压制后的信号，取得一定效果。

EMD方法不需要设置针对性的参数，具有较强的泛化能力。但如果在频率—空间域的频率切片体上进行EMD处理，边界效应等现象常常导致分解过程精度低或不稳定^[25]。在EMD方法的基础上发展的变分模态分解(VMD)方法，通过将地震数据进行解析化，在时频域内将地震信号分解为各个模态分量。对于时频域中带限频率切片的处理，VMD方法^[26]明显优于传统EMD方法，是一种自适应、完全非递归的模态变分方法。基于奇异值分解的VMD方法通过确定模态分解个数，自适应地匹配每种模态的最优中心频率和有限带宽，最终得到变分问题的最优解，实现IMF的高效分离。该方法克服了传统EMD的端点效应与模态分量混叠问题，求解过程具有更高的稳定性。

本文通过希尔伯特-黄变换(HHT)将地震数据解析化，在时频域进行变分模态分解。通过分频分解的思路，分析每个模态分量各频段上的噪声和有效信号的能量分布，设置合适的分量数，分频段压制噪声，再将有效信号模态分量重构回时频域。实际资料测试结果表明，相较于常规RNA随机噪声压制方法，该方法有效地压制了随机背景噪声，同时对陡倾角的线性干扰有明显的压制效果。

1 算法原理

对于地震剖面中的线性同相轴，可将其描述为空间坐标x和时间t的函数

$ u(t, x)=w\left(t-\frac{x}{c}\right) $

(1)

式中：w为地震子波；c为波的视速度。在频率域可由Fourier变换表示为

$ \hat{u}(f, x)=W(f) \exp \left(\frac{\mathrm{i} 2 \pi f x}{c}\right) $

(2)

式中W(f)为子波的频谱。

上式表明$\hat{u}(f, x) $可视为由一组复谐波叠合而成的信号。在此基础上，可将目标数据的解析信号直接进行分解，由此构造的优化问题目标函数及其增广拉格朗日函数分别为

$ \begin{array}{*{20}{c}} &\mathop {\min }\limits_{\left\{ {{u_k}} \right\}\langle {\omega _k}\} } \sum\limits_{k=1}^{K}\left\|\partial x\left\{\left[\delta(x)+\frac{\mathrm{i}}{\pi x}\right] u_{k}(x) \exp \left(-\mathrm{j} \omega_{k} x\right)\right\}\right\|_{2}^{2}\\ &\text { 满足 } \;\;\; \sum\limits_{k=1}^{K} u_{k}=u \end{array} $

(3)

$ \begin{array}{l} L\left(\left\{u_{k}\right\},\left\{\omega_{k}\right\}, \lambda\right)=\alpha \sum\limits_{k=1}^{K} \| \partial_{x}\left[u_{k}(x)\right] \times \\ \;\;\;\;\; \exp \left(-\mathrm{j} \omega_{k} x\right)\left\|_{2}^{2}+\right\| u(x)-\sum\limits_{k=1}^{K} u_{k}(x) \|_{2}^{2}+ \\ \;\;\;\;\; \left\langle\lambda(x), u(x)-\sum\limits_{k=1}^{K} u_{k}(x)\right\rangle \end{array} $

(4)

式中：α为权重系数，调节目标函数中的变分项与保真项；λ(x)为拉格朗日算子；$\left[\delta(x)+\frac{\mathrm{j}}{\pi x}\right] u_{k}(x) $为原始数据u(x)对应的解析信号；ω_k为k阶模态分量的中心频率；K为模态分量总阶数；〈·，·〉表示两个信号的内积。对于目标函数的最优化求解，较为理想的方法可采用交替方向乘子法(ADMM)^[27]，采用线性逼近的思路，通过调节线性项系数逐渐接近最优解，使得最优化问题在二次惩罚项系数很小的情况下，仍可得到满足精度要求的解。其核心思路是对中心频率和相应模态的更新

$ \begin{aligned} u_{k}^{n+1}=& \arg \min _{u_{k}} L\left(\left\{u_{i<k}^{n+1}\right\},\left\{u_{i \geqslant k}^{n}\right\},\left\{\omega_{i}^{n}\right\}, \lambda^{n}\right) \\ =& \arg \min _{u_{k}} \alpha\left\|\partial_{x}\left[u_{k}(x) \mathrm{e}^{-\mathrm{j} \omega_{k} x}\right]\right\|_{2}^{2}+\\ &\left\|(x)-\sum\limits_{i} u_{i}(x)+\frac{\lambda(x)}{2}\right\|_{2}^{2} \end{aligned} $

(5)

对上式进行傅里叶变换后，可得频率域解

$ {\hat u_k^{n + 1}(\omega ) = \frac{{\hat u(\omega ) - \sum\limits_{i \ne k} {{{\hat u}_i}} (\omega ) + \frac{{\hat \lambda (\omega )}}{2}}}{{1 + \frac{\alpha }{2}{{\left( {\omega - {\omega _k}} \right)}^2}}}} $

(6)

$ {\omega _k^{n + 1} = \arg \mathop {\min }\limits_{{\omega _k}} L\left( {\left\{ {u_i^{n + 1}} \right\},\left\{ {\omega _{i < k}^n} \right\},\left\{ {\omega _{ik}^n} \right\},{\lambda ^n}} \right)}\\ \;\;\;\; { = \frac{{\int_{ - \infty }^\infty \omega {{\left\| {{{\hat u}_k}(\omega )} \right\|}^2}{\rm{d}}\omega }}{{\int_{ - \infty }^\infty {{{\left\| {{{\hat u}_k}(\omega )} \right\|}^2}} {\rm{d}}\omega }}} $

(7)

式中：$\hat{u}_{k} $为u_k的傅里叶变换；n为更新次数。

算法基本流程如下：

(1) 应用合适时窗选取部分目标数据，进行VMD参数初始化；

(2) 应用HHT将目标数据从时间域变换到时频域；

(3) 对各个频率切片进行VMD分解，确定各模态分量的中心频率(图 1)；

(4) 按能量强弱对VMD分解后得各模态分量(记为VMDIMF)排序，保留前若干个有效VMD-IMF，重构得到滤波后的频率切片；

(5) 将滤波后的所有频率切片反变换回时空域；

(6) 移动时窗，重复上述步骤，直到完成所有数据的处理。

对于三维数据体，仅增加一项正交空间坐标基即可，整体算法流程不变。

对地震数据进行HHT，可分析含噪非平稳信号在不同频率段内能量沿时间轴的分布。图 2为原始单道地震记录及其前半段(1~500ms)时频谱，其中IMF1~IMF6依次对应信号的从高频至低频信息(对于单道数据，频率切片简化为IMF)。由图 2可见，高频端和低频端主要由噪声构成，有效信息集中在IMF3~IMF5。以IMF5的时频谱为例，对其进行VMD处理，分解得到各阶VMDIMF，如图 3所示。图中将各VMDIMF的能量进行归一化，方便对比分析。若分量中存在异常值，可对图 2进行综合分析，决定是将其作为异常噪声进行剔除，还是作为局部有效信号予以保留。

用2、25、288Hz谐波合成信号对比VMD法与传统EMD法分解效果(图 4和图 5)。在分解效率方面，由于VMD递归确定最优中心频率，仅3个VMD-IMF即分解出三个频率成分，而传统EMD方法需要6个IMF才分解出288Hz成分；在精度方面，VMD法分解得到的VMDIMF2、VMDIMF3、VMDIMF4均准确定位了合成信号的各频率成分(图 4)，而EMD法分解结果中IMF2、IMF3与IMF5均为无效分量(图 5)。可见，对谐波合成信号的分解，VMD法具有明显的效率和精度优势。

2 模型测试

使用VMD法对地震资料进行处理时，首先需要指定总的模态个数K，该参数直接影响去噪效果。在实际应用前，首先对小块资料(比如200×200个样点)进行测试。对模型数据(图 6a)加入频带为8~55Hz的高斯随机噪声，噪声平均能量大于有效信号能量的50%，其信噪比为3.09dB(图 6b)。该加噪数据不同K值的VMD法去噪结果如图 7所示。当K=3时背景噪声压制效果最好，但有效成分出现较为明显的割裂痕迹(图 7a红色箭头所示)，是由时窗截取的边缘效应被过低K值放大所致，可根据实际情况，调整时窗长度。在K为5或7时，去噪结果的割裂痕迹明显得到改善；随着K值增大，更多的背景噪声被保留，当K值过大时(该例K大于7时)，出现纵向高倾角处理痕迹(图 7d~图 7f黄框所示)。图 8为图 7去噪后数据的时频谱，图 9为去噪结果信噪比(SNR)随K值的变化曲线，可见：K值较小时，处理结果的信噪比更高；随着K值的增大，更多噪声成分被保留，信噪比逐渐下降。

3 实际资料应用

实际地震剖面如图 10a所示，采样间隔为2ms，记录长度为8s，道间距为12.5m，具有较强的随机噪声干扰，还有部分线性相关噪声，中深层同相轴连续性较差。本文方法去噪结果如图 10b所示，与原始剖面相比，浅层水平同相轴和深层低伏度同相轴的质量得到明显提升，连续性得到明显加强，构造形态更加清晰。与常规RNA去噪方法的处理结果(图 10c)相比，本文方法去噪结果未出现轻微类似斜干扰的人为痕迹(红色椭圆所示)。通过局部放大对比(图 11、图 12)，处理前浅层水平同相轴背景有残留的面波成分(图 11a)经处理后得到明显的压制(图 11b)；中深层的低幅度构造，在处理前受到线性噪声和高频背景噪声干扰(图 12a)，处理后线性噪声痕迹几乎消失，背景随机噪声能量得到较强压制，与有效同相轴的能量不在同一数量级(图 12b)。由处理前、后的频谱(图 13)可见，由于压制了低频面波和高频背景噪声能量，低频能量和高频能量有所减弱。

4 结束语

本文通过HHT得到地震资料在各频率切片上的能量分布，然后通过VMD对各频率切片进行分解；分析有效信号与噪声在时频切片上的能量分布，在此基础上优选出有效信号模态分量重构时频切片；最后反变换回时空域，达到噪声压制的目的。分频处理可以在不同频率段内对二维信号，甚至单道信号进行有效信息的筛选，从而对局部异常噪声进行更精细的压制。模型数据测试和实际资料应用结果表明，与常规RNA方法相比，本文方法可有效压制线性干扰和随机背景噪声，处理后资料信噪比更高。