2020, Vol. 63
地震勘探中的有效信号总是不可避免地被随机噪声所干扰, 尤其在陆上勘探中, 随机噪声更加严重地影响有效信号的准确成像.一方面, 随着国家深地资源勘查开采计划的开展, 对深层地震数据的质量有了更高的需求, 深层地震数据的低信噪比问题是随机噪声和有效信号能量衰减共同作用的结果.另一方面, 当前油气勘探的目标已经转变为“两宽一高”, 更大规模的数据量给随机噪声压制方法提出了更加苛刻的处理效率问题.另外, 人工地震数据本质上是时空变化的, 表现在地震同相轴的能量、轨迹、时频谱等特征随时间和空间位置不同而变化(张雅晨等, 2019), 传统的随机噪声处理方法很难适应实际有效信号的时变特征.因此, 开发快速、有效的随机噪声衰减方法仍然是当前人工地震勘探数据处理亟需解决的核心问题之一.当前, 地震随机噪声衰减方法的不同研究方向包括：基于随机噪声统计学特征的均值滤波类和中值滤波类方法(Bonar and Sacchi, 2012; Liu et al., 2009; Liu, 2013; Zhao et al., 2016), 基于各种数学变换的随机噪声压制方法, 如seislet变换、bandelet变换、双树复小波域双变量方法(刘洋等, 2009; Wang et al., 2015; 汪金菊等, 2016)等, 基于数学分析工具的处理方法, 如经验模态分解、奇异值分解(Bekara and Van Der Baan, 2009; 刘志鹏等, 2012), 基于数据轨迹特征的处理方法, 如边缘保护滤波、方向可控滤波(Lu and Lu, 2009; 黄梅红和李月, 2016), 基于反演的去噪方法, 如贝叶斯反演, 基于反演的时空域去噪方法(Yuan et al., 2012; Zhao et al., 2017)等, 以及基于预测滤波类的随机噪声处理方法, 如自回归滑动平均滤波(ARMA)、正则化非平稳自回归滤波(RNA)、流式正交预测滤波(Sacchi and Naghizadeh, 2009; Liu et al., 2015; Liu and Li, 2018)等.在不同类型的随机噪声衰减方法中, 预测滤波方法仍然是工业界应用的主流技术.
地震同相轴在空间方向具有可预测性, 这种预测特征能够用于分离强随机噪声中的目标同相轴, 该类方法可以在频率-空间域(f-x域)或时间-空间域(t-x域)内实现.f-x域预测滤波技术由Canales(1984)提出并且由Gulunay(1986)进一步发展, 已经成为一种工业界标准的方法——FXDECON.Claerbout(1992)提出预测滤波过程也可以在t-x域中实现.Abma和Claerbout(1995)讨论了在f-x域和t-x域预测线性同相轴方法的不同, FXDECON滤波器相当于将时间方向上的滤波器长度选取为整道的t-x域预测滤波器, FXDECON相比较t-x域方法在去噪能力方面略差, 但是具有更快的计算速度, 这在三维地震数据处理中具有更大的优势.常规的FXDECON方法只能处理倾角线性变化的同相轴信息, 为了处理空间方向弯曲的非平稳同相轴, Liu等(2012)提出了f-x域RNA方法用于压制地震随机噪声, 该方法通过正则化条件使得非平稳滤波器系数在空间和频率方向平滑, 并且将其拓展到三维情况(Liu et al., 2013).尽管基于RNA的f-x域非平稳预测滤波方法可以提供很好的自适应滤波特征, 但是该方法使用迭代递归算法来计算目标函数中随频率-空间变化的滤波器系数, 这导致了计算成本增加, 尤其在三维宽方位数据处理中, 变系数的存储极大地限制这种方法的广泛应用.Fomel和Claerbout(2016)提出了t-x域流式预测误差滤波器(PEF)的概念, 该方法在新数据通过滤波器时更新PEF的系数, 将滤波器系数的计算成本降低到单个卷积的数量级, 并用于解决缺失数据的插值问题.Liu和Li(2018)进一步提出了t-x域流式正交预测滤波器, 为自适应预测滤波压制随机噪声提供了快速的求解思路.
本文基于流式处理框架提出一种新的f-x域流式预测滤波方法, 通过建立频率域的预测自回归方程, 利用谢尔曼-莫里森(Sherman-Morrison)公式直接推导得到最小二乘目标函数在频率域中的非迭代解, 利用求得的预测滤波系数可以恢复随机噪声环境中的局部平面波信息.该方法在同一频率切片下, 利用空间方向上相邻滤波器的相似性实现数据的预测, 滤波器虽然具有固定长度, 但滤波器系数在预测过程中被不断更新.理论模型和实际数据的测试结果说明, 相比工业标准的FXDECON方法和f-x域RNA方法, 本文提出的f-x域流式预测滤波技术能够在计算效率和处理效果方面达到更好的平衡.
1 f-x域预测滤波基础理论 1.1 FXDECON预测滤波理论假设地震数据s是关于时间t和空间位置x的平面波函数.地震数据可以表示为N个斜率为pi的倾斜同相轴的线性叠加, 即：
|
(1) |
其中, wi是与第i个反射同相轴相关的地震子波, 与位于时间ti的狄拉克函数进行褶积, 并且ti=τi+pix, τi是某个同相轴在x=0位置处的截距时间, 将(1)式进行傅里叶变换可得
|
(2) |
Wi(ω)是子波wi(t)的傅里叶变换.根据自回归(AR)模型, 频率域第k个地震道S(ω, xk)可以表示为不同空间位置地震道S(ω, xk－i)的叠加：
|
(3) |
其中, 
公式(3)可以写成f-x域自回归方程：
|
(4) |
每一个频率切片中, 未知数的个数为N, 公式(4)可以写成XC=S的复数矩阵形式, 通常用迭代的方法求解滤波器系数, 如FXDECON方法, 将预测误差表示为：r=S－XC, 误差能量为：I=(S－XC)T(S－XC), 其中T表示复数共轭转置.使其导数等于零将误差能量最小化：
|
(5) |
将(5)式表示为：
|
(6) |
定义自相关矩阵R=XTX和互相关矩阵g=XTS.其中R=P+jQ, P是实数域上的对称矩阵, Q是实数域上的斜对称矩阵, 将(6)式写为如下形式：
|
(7) |
其中Re和Im分别代表实部和虚部, 可以通过迭代的高斯消元法计算滤波器系数C.
1.2 f-x域流式预测滤波理论如果地震同相轴为非线性形态, 可以假设在空间方向存在局部线性关系, 此时为了能够预测非平稳变化的地震同相轴, 公式(3)可以修改为：
|
(8) |
其完整的矩阵形式可以写为：
|
(9) |
其中, 任意数据点的自回归方程可以表达为：
|
(10) |
上述数学欠定问题可以通过增加正则化条件进行求解(Liu et al., 2013), 此时每个数据样点的预测滤波因子依赖于空间坐标.Fomel和Claerbout(2016)提出一种流式处理算法代替正则化条件, 本文方法基于这种流式框架, 当计算新的数据点S(ω, xk)的滤波器系数时, 假设新的预测滤波器与相邻的滤波器在空间方向上是相似的, 通过增加一个能控制滤波器变化的限制条件, 使得新的滤波器系数Ci(ω, xk)逼近于先前的滤波器系数
|
(11) |
公式(11)可以表示为：
|
(12) |
其中, I代表单位矩阵, S(ω, x)是在f-x域内(ω, x)点位置的数据, λ为用于调节C偏离
|
(13) |
方程(12)在复数域的最小二乘解为：
|
(14) |
由于矩阵逆的求解需要大量的运算, Fomel和Claerbout(2016)使用线性代数中的谢尔曼-莫里森(Sherman-Morrison)公式直接求解逆矩阵问题以避免迭代计算.该公式是一种解决特殊矩阵逆的解析算法(Hager, 1989).如果矩阵A, I－VA-1U和A-UV都是可逆的, A-UV的逆矩阵可以由以下公式计算：
|
(15) |
在特殊的情况下, 当矩阵U是列向量u, V是行向量v, 方程(15)可以被写为：
|
(16) |
利用谢尔曼-莫里森公式可以将公式(14)中的逆矩阵转化为：
|
(17) |
得到矩阵的最小二乘解为：
|
(18) |
化简得到滤波器系数为：
|
(19) |
公式(19)表明可以通过加上按照比例缩放的数据来自适应更新滤波器系数, 该过程仅需要简单的复数运算.最后, 可预测的信号和不可预测的随机噪声分别表示为：
|
(20) |
|
(21) |
将公式(19)代入公式(20)和(21)中, 化简得到：
|
(22) |
|
(23) |
本方法假设随机噪声在频率-空间域仍然表现为随机特征, 计算中仅需要存储前一个空间样点上的滤波器系数
为了验证本文方法的有效性, 选取标准叠后地震数据模型, 模型中包含上层倾斜同相轴, 下层以正弦规律变化的弯曲同相轴, 一个不整合面和一个倾斜的断层(图 1a), 图 1b为加入了高斯随机噪声后的模型数据.为了验证所提出方法的有效性, 我们选取工业标准的FXDECON方法与其对比.图 1c为用FXDECON处理的结果(信噪比为8.493), 可以看出滤波后的剖面中仍然含有少量随机噪声, 差剖面(图 1d)显示该方法压制随机噪声的同时也去除了较多有效信号, 对构造信息的保护能力较差.为了提高处理效果, 往往将数据分割为小的处理窗口, 以满足在每个子窗口内同相轴为线性的假设条件, 但是这种方法面临窗口大小的自适应选择难题.f-x域RNA方法滤波结果如图 1e所示(信噪比为11.62), 滤波效果较好, 大量随机噪声被压制, 从差剖面(图 1f)中看出该方法对弯曲同相轴的构造信息保护较好, 仅滤除少量有效信号.本文提出的f-x域流式预测滤波处理结果如图 1g所示(信噪比为9.747), 通过计算得到Smax值为0.129, 选取滤波器参数λ为0.13, 滤波器长度为25.从图中可以看到同相轴连续性增强, 剖面整体清晰度增加, 但有少量微弱的虚假同相轴存在, 这种现象在所有的f-x域预测滤波方法中都存在, 可以通过增加频率方向的滤波器系数平滑约束得以改善.图 1h为f-x域流式预测滤波器分离的随机噪声, 通过图 1d与图 1h的对比, 本文方法能够在有效压制随机噪声的同时更好地保护时空变信号.通过图 1f与1h对比, 本文方法对倾斜同相轴信息保护效果更好, 但对弯曲构造保护效果略差.这是由于f-x域RNA方法能更好地预测出局部变化特征, 而f-x域流式预测滤波对全局特征预测效果更好.同时, 本文方法的计算效率要优于基于迭代的自适应预测滤波方法.
|
图 1 理论模型测试结果 (a)理论模型数据; (b)加入随机噪声后数据; (c) FXDECON去噪结果; (d) FXDECON分离的噪声; (e) f-x域RNA去噪结果; (f) f-x域RNA分离的噪声; (g) f-x域流式预测滤波去噪结果; (h) f-x域流式预测滤波分离的噪声. Fig. 1 Synthetic data test (a) Model; (b) Noisy data; (c) The denoised result by FXDECON; (d) The noise removed by FXDECON; (e) The denoised result by f-x RNA; (f) The noise removed by f-x RNA; (g) The denoised result by f-x streaming prediction filter; (h) The noise removed by f-x streaming prediction filter. |
计算图 1中理论模型及处理结果的f-k谱, 进一步验证本文方法的有效性.图 2a为原始信号的f-k谱, 与理论模型有很好的对应性.两侧能量较强的信号对应着深层以正弦形式变化的弯曲同相轴, 右侧斜率较大的能量较强的信号为浅层的倾斜同相轴.图 2b为加入高斯随机噪声后模型数据的f-k谱, 从图中看出随机噪声在f-k域中分布也是随机的.FXDECON方法、f-x域RNA方法和f-x域流式预测滤波方法处理后数据的f-k谱分别如图 2c, 2e和2g所示.可以看到本文方法滤波效果合理, 去除了大量的噪声, 滤波后谱特征与原始理论模型相似度更高.图 2d, 2f和2h分别为使用FXDECON方法、f-x域RNA方法和f-x域流式预测滤波方法滤除的随机噪声部分, 可以看出所提出方法对有效信号的保护能力和压制随机噪声能力优于传统FXDECON方法, 对倾斜同相轴构造的保护效果优于f-x域RNA方法.
|
图 2 理论模型f-k谱 (a)原始信号f-k谱; (b)含噪数据f-k谱; (c) FXDECON滤波剖面f-k谱; (d) FXDECON压制噪声f-k谱; (e) f-x域RNA滤波剖面f-k谱; (f) f-x域RNA压制噪声f-k谱; (g) f-x域流式预测滤波剖面f-k谱; (h) f-x域流式预测滤波压制噪声f-k谱. Fig. 2 f-k spectra of synthetic data (a) f-k spectra of model; (b) f-k spectra of noisy data; (c) f-k spectra of the denoised result by FXDECON; (d) f-k spectra of the noise removed by FXDECON; (e) f-k spectra of the denoised result by f-x RNA; (f) f-k spectra of the noise removed by f-x RNA; (g) f-k spectra of the denoised result by f-x streaming prediction filter; (h) f-k spectra of the noise removed by f-x streaming prediction filter. |
在实际数据测试中, 我们选取某地区时间偏移剖面, 该数据共310个地震道, 在时间方向上有700个采样点, 采样间隔为1 ms(图 3).数据中存在比较复杂的构造信息, 包括浅层的平行同相轴和深层的倾斜同相轴.使用FXDECON与f-x域流式预测滤波进行对比, FXDECON的处理结果如图 4a所示, 从图中看出部分随机噪声被压制, 剖面整体清晰度增加, 浅层同相轴连续性增强.将滤波后的结果与原始数据作差得到滤除的噪声, 如图 4b所示.在噪声剖面中可以看到标准FXDECON方法对非平稳有效信号振幅的保护效果差, 有效信号损失较为严重, 同时随机噪声压制效果有限.接下来, 选取f-x域RNA的方法得到如图 4c的滤波结果, 可以看出大量随机噪声被压制, 从差剖面(图 4d)中看出该方法对构造信息的保护效果较好, 仅滤除少量有效信号.通过计算得到Smax值为297.3, 选取参数λ为350, 空间滤波长度为30的f-x域流式预测滤波器对实际数据进行处理(图 4e), 可以看出地震同相轴的连续性和局部构造特征得到提高, 剖面整体清晰度增加, 同FXDECON相比, 噪声剖面(图 4f)中含有更少的有效信号.同f-x域RNA方法相比, 虽然去噪效果略差, 但计算效率显著提高.另外, 可以通过调整参数λ的值在随机噪声压制和有效信号保护之间进行取舍.该方法也可以对三维叠前地震数据进行处理, 即对每一个单炮数据以二维形式处理, 或将本文方法拓展为空间二维形式, 即在公式(12)中增加沿空间y方向的约束项.
|
图 3 实际叠后数据 Fig. 3 Poststack field data |
|
图 4 实际数据测试结果 (a) FXDECON去噪结果; (b) FXDECON去除的噪声; (c) f-x域RNA去噪结果; (d) f-x域RNA去除的噪声; (e) f-x域流式预测滤波去噪结果; (f) f-x域流式预测滤波去除的噪声. Fig. 4 Field data test (a) The denoised result by FXDECON; (b) The noise removed by FXDECON; (c) The denoised result by f-x RNA; (d) The noise removed by f-x RNA; (e) The denoised result by f-x streaming prediction filter; (f) The noise removed by f-x streaming prediction filter. |
接下来, 计算野外实际数据的f-k谱来验证本文方法的有效性.图 5a为实际数据的f-k谱, 随机噪声的频率集中在150 Hz以内, 在f-k域中随机分布.FXDECON方法处理后数据的f-k谱如图 5b所示, 可以看出仍有部分随机噪声没有被滤除, 损失了部分有效信号, 对构造信息保护效果较差.f-x域RNA方法处理后数据的f-k谱如图 5c所示, 可以看出随机噪声被有效地压制.f-x域流式预测滤波方法处理后数据的f-k谱如图 5d所示, 与图 5b和图 5c的对比可以看出本文方法压制随机噪声的能力优于FXDECON方法, 但与f-x域RNA方法相比效果略差.
|
图 5 实际叠后数据f-k谱 (a)实际叠后数据f-k谱; (b) FXDECON滤波剖面f-k谱; (c) f-x域RNA滤波剖面f-k谱; (d) f-x域流式预测滤波剖面f-k谱. Fig. 5 f-k spectra of poststack field data (a) f-k spectra of poststack field data; (b) f-k spectra of the denoised result by FXDECON; (c) f-k spectra of the denoised result by f-x RNA; (d) f-k spectra of the denoised result by f-x streaming prediction filter. |
本文基于流式处理框架提出一种新的f-x域流式预测滤波方法, 建立频率域的预测滤波自回归方程, 利用复数逆矩阵直接求解复数域最小二乘问题, 非迭代计算时空变化的滤波器系数.该滤波方法构建流动式计算过程不断更新滤波器系数, 有效降低处理非平稳地震数据的复杂度.通过与传统方法如工业标准FXDECON方法和适用于处理非平稳地震数据的迭代f-x域RNA方法相比, 本文方法用近似解析解代替传统迭代最小二乘解有效地提高运算速度, 能更好地平衡有效信号保护、随机噪声压制以及计算效率的关系, 得到合理的滤波结果.理论模型和实际数据的测试结果表明提出的方法能够在降低计算成本的同时有效地衰减非平稳地震数据中的随机噪声.
致谢 感谢中国石油大学(北京)刘国昌副教授的讨论, 感谢三位匿名审稿专家提出的宝贵意见.
Abma R, Claerbout J. 1995. Lateral prediction for noise attenuation by T-X and F-X techniques. Geophysics, 60(6): 1187-1896. |
Bekara M, Van Der Baan M. 2009. Random and coherent noise attenuation by empirical mode decomposition. Geophysics, 74(5): V89-V98. DOI:10.1190/1.3157244 |
Bonar D, Sacchi M. 2012. Denoising seismic data using the nonlocal means algorithm. Geophysics, 77(1): A5-A8. DOI:10.1190/geo2011-0235.1 |
Canales L L. 1984. Random noise reduction.//54th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 525-527.
|
Claerbout J F. 1992. Earth Soundings Analysis:Processing Versus Inversion. Oxford: Blackwell Scientific Publications.
|
Fomel S, Claerbout J. 2016. Streaming prediction-error filters.//86th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 4787-4791.
|
Gulunay N. 1986. Fxdecon and complex wiener prediction filter.//56th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 279-281.
|
Hager W W. 1989. Updating the inverse of a matrix. SIAM Review, 31(2): 221-239. DOI:10.1137/1031049 |
Huang M H, Li Y. 2016. Suppression of seismic random noise based on steerable filter. Chinese Journal of Geophysics (in Chinese), 59(5): 1815-1823. DOI:10.6038/cjg20160524 |
Liu G C, Chen X H, Du J, et al. 2012. Random noise attenuation using f-x regularized nonstationary autoregression. Geophysics, 77(2): V61-V69. DOI:10.1190/geo2011-0117.1 |
Liu G C, Chen X H. 2013. Noncausal f-x-y regularized nonstationary prediction filtering for random noise attenuation on 3D seismic data. Geophysics, 93: 60-66. |
Liu Y, Fomel S, Liu C, et al. 2009. High-order seislet transform and its application of random noise attenuation. Chinese Journal of Geophysics (in Chinese), 52(8): 2142-2151. DOI:10.3969/j.issn.0001-5733.2009.08.024 |
Liu Y, Liu C, Wang D. 2009. A 1D time-varying median filter for seismic random, spike-like noise elimination. Geophysics, 74(1): V17-V24. DOI:10.1190/1.3043446 |
Liu Y, Liu N, Liu C. 2015. Adaptive prediction filtering in t-x-y domain for random noise attenuation using regularized nonstationary autoregression. Geophysics, 80(1): V13-V21. |
Liu Y, Li B X. 2018. Streaming orthogonal prediction filter in t-x domain for random noise attenuation. Geophysics, 83(4): F41-F48. DOI:10.1190/geo2017-0322.1 |
Liu Y K. 2013. Noise reduction by vector median filtering. Geophysics, 78(3): V79-V87. DOI:10.1190/geo2012-0232.1 |
Liu Z P, Zhao W, Chen X H, et al. 2012. Local SVD for random noise suppression of seismic data in frequency domain. Oil Geophysical Prospecting (in Chinese), 47(2): 202-206. |
Lu Y H, Lu W K. 2009. Edge-preserving polynomial fitting method to suppress random seismic noise. Geophysics, 74(4): V69-V73. DOI:10.1190/1.3129907 |
Sacchi M D, Naghizadeh M. 2009. Adaptive linear prediction filtering for random noise attenuation.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 3347-3351.
|
Wang J J, Yuan L, Liu W R, et al. 2016. Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation. Chinese Journal of Geophysics (in Chinese), 59(8): 3046-3055. DOI:10.6038/cjg20160827 |
Wang X K, Gao J H, Chen W C, et al. 2015. The seismic random noise attenuation method based on enhanced bandelet transform. Journal of Applied Geophysics, 116: 146-155. DOI:10.1016/j.jappgeo.2015.03.002 |
Yuan S Y, Wang S X, Li G F. 2012. Random noise reduction using Bayesian inversion. Journal of Geophysics and Engineering, 9(1): 60-68. DOI:10.1088/1742-2132/9/1/007 |
Zhang Y C, Liu Y, Liu C, et al. 2019. Seismic random noise attenuation using FDOC-seislet transform and threshold for seismic data with low SNR. Chinese Journal of Geophysics (in Chinese), 62(3): 1181-1192. DOI:10.6038/cjg2019M0062 |
Zhao X, Li Y, Zhuang G H, et al. 2016. 2-D TFPF based on Contourlet transform for seismic random noise attenuation. Journal of Applied Geophysics, 129: 158-166. DOI:10.1016/j.jappgeo.2016.03.030 |
Zhao Y M, Li G F, Wang W, et al. 2017. Inversion-based data-driven time-space domain random noise attenuation method. Applied Geophysics, 14(4): 543-550. DOI:10.1007/s11770-017-0647-4 |
黄梅红, 李月. 2016. 基于方向可控滤波的地震勘探随机噪声压制. 地球物理学报, 59(5): 1815-1823. DOI:10.6038/cjg20160524 |
刘洋, Fomel S, 刘财, 等. 2009. 高阶seislet变换及其在随机噪声消除中的应用. 地球物理学报, 52(8): 2142-2151. DOI:10.3969/j.issn.0001-5733.2009.08.024 |
刘志鹏, 赵伟, 陈小宏, 等. 2012. 局部频率域SVD压制随机噪声方法. 石油地球物理勘探, 47(2): 202-206. |
汪金菊, 袁力, 刘婉如, 等. 2016. 地震信号随机噪声压制的双树复小波域双变量方法. 地球物理学报, 59(8): 3046-3055. DOI:10.6038/cjg20160827 |
张雅晨, 刘洋, 刘财, 等. 2019. 低信噪比地震资料FDOC-seislet变换阈值消噪方法研究. 地球物理学报, 62(3): 1181-1192. DOI:10.6038/cjg2019M0062 |