2017, Vol. 60
信噪比的高低是衡量地震资料质量的重要指标之一，也是评价去噪方法优劣的重要指标.如何准确估算地震信号的信噪比一直是有待解决的地震数据处理问题.在信噪比估值方法中，信噪比计算公式的选取起着重要作用，一个稳健的计算公式能得到可信的估值结果.近几十年来，研究人员提出了多种不同的信噪比计算公式.Turin (1960)将信号相关峰值期望平方和峰值方差的比值做为光学图样中匹配滤波器的信噪比定义.Kerr (1966)提出了把信号功率和高斯白噪方差的平方做商的最大似然估计法.Pratt (1978)提出了利用信号图像的边缘亮度落差和噪声方差来估计信噪比的计算公式.另外，因为信号的各阶统计量代表了信号不同的统计学特征，所以国外还涌现出一批以信号的多阶统计量来定义信号信噪比的研究学者.Benedict和Soong (1967)首先使用一阶统计量和二阶统计量的比值求信噪比.Matzner (1993)把该方法提升到二阶四阶统计量的计算.Loópez-Valcarce和Mosquera (2007)认为传统的二四阶统计量方法不适用高信噪比的计算，因此提出一种六阶统计量计算方法，该方法通过使信噪比方差最小化来添加约束条件，然后迭代求解信噪比.Álvarez-Díaz等 (2010)则提出一种综合利用二到八偶数阶的信噪比计算公式并给出两种最优化求解方法，以此进一步提高计算的准确性和有效适用范围.上述信噪比计算公式在各自领域都得到了很好的计算结果，然而它们在应用于地震信号时会出现信号特性不匹配或因为失去特定的约束条件而无法求解的问题.近年来，地震信号的信噪比公式没有得到太多研究，目前常用的是信号全局能量与噪声全局能量的比值.因此如何提出一个适用于非平稳地震信号的更稳健信噪比公式是亟需解决的问题.
在信噪比计算公式中，估计出有效信号是计算的核心.目前已有的地震信号信噪比估值方法在有效信号估计上取得了较多的研究成果.霍志坚和傅旦丹 (1990)提出利用地震信号各道叠加求均值和最大相关峰值时移量来提取信号的统计平均分析法.刘洋和李承楚 (1997)提出了叠加法和SVD法两种信号估计方法.赵亮等 (2003)利用信号互相关而噪声互不相关的原理来提取有效信号功率谱.牛聪等 (2006)利用信号和噪音对应频率成分不同的原理提出了频率域信噪比估值法，同时提出了考虑信号相位变化的相关时移估值法.王红玲 (2007)在相关时移法的基础上利用F-K滤波分离相关干扰，以此提高地震信号相关时移的准确度.张军华等 (2009)综合比较了几种常用估算方法的计算效果并给出了每种方法的适用情况.
本文提出了一种局部信噪比计算公式，主要用于表征非平稳地震数据中信号与噪声的分布关系.其基本原理是将信噪比计算问题转化为反演问题，对瞬时信噪比计算使用正则化共轭梯度法迭代寻找最优解，然后通过正则化参数的选取来求得信号中不同位置数据点的局部信噪比.本文使用Tikhonov正则化 (Tikhonov and Arsenin, 1977) 条件计算数据平滑度固定不变的局部正则化算子，求出具有相同平滑特性的局部信噪比.理论模型与实际数据验证的结果表明，相比传统的信噪比估值方法，局部信噪比估值法能更直观准确地反映地震数据非平稳变化区域中信噪比的大小和特征.
2 理论基础 2.1 局部信噪比公式假设地震数据可表示为信号与噪声的叠加：
|
(1) |
其中，d是待求信噪比数据, s是d的真实有效信号，n是噪声.传统信噪比的定义是信号和噪声某种属性的比值，常规地震信号全局信噪比计算公式可以表示为
|
(2) |
传统的信噪比计算方法是将地震信号中所有数据振幅求平方和以计算全局信噪比，这种方法只能评价整体地震资料质量好坏，无法给出不同数据点位置的信噪比信息.因此，可以考虑在地震信号中截取时空窗口，得到窗口化形式的信噪比计算公式：
|
(3) |
其中，m是窗口长度.SNR (t, x) 是时间为t，空间为x的数据点处的窗口信噪比.当m>1时，公式 (3) 表示以时窗中所有数据的平方和来近似点 (t, x) 处的信噪比，此时窗口长度越大，计算出的信噪比越不准确.当m=1时，时窗中只有一个数据点，公式 (3) 表示地震数据的点点相除，此时计算出的信噪比在理论上来讲是最准确的，本文称之为“瞬时信噪比”.实际处理中含随机噪声的地震信号存在很多零值点，因此把数据直接相除是不可取的.可以在计算中使有效信号为零的数据点信噪比等于零，在对数域中令其等于全局最低信噪比；对有效信号不为零而噪声为零的点，赋给噪声一个小的扰动量使得除法有解，得到点点相除的信噪比结果.然而这种近似方法本身计算的结果不准确，再加上地震信号中异常值点较多，最后求出的局部信噪比剖面往往没有参考意义.本文提出一种局部信噪比计算公式来代替公式 (3) 所示的窗口化信噪比计算公式.
定义瞬时信噪比计算公式如下：
|
(4) |
其中，
|
(5) |
n(i, j)=d(i, j)－s(i, j).可以利用最小二乘原理构建如下目标函数求解公式 (5)：
|
(6) |
公式 (6) 对应的最小二乘解为
|
(7) |
其中，N=diag[n(i, j)2]，是以n(i, j)2为主对角线元素的矩阵.S=diag[s(i, j)2].如果矩阵NTN可逆，则R可解.如果矩阵NTN不满秩，可以加入正则化条件，构建如下目标函数：
|
(8) |
其中，D是正则化算子.ε是惩罚因子，用来平衡最小二乘目标函数和正则化条件的权重.正则化条件下的最小二乘解为
|
(9) |
正则化算子D控制着目标函数的局部平滑特性，不同的正则化算子得到不同的平滑结果.本文选取Tikhonov正则化作为正则化算子，正则化算子的长度可以与m相等.Tikhonov正则化是对数据域的正则化，可以保证拟求取的信噪比具有一致的平滑参数，使计算的局部信噪比具有一致的衡量标准.
公式 (9) 可以通过共轭梯度算法来寻找最优解，该最优解为局部平滑解，最终局部正则化信噪比的计算公式为
|
(10) |
公式 (10) 可以写为如下只含有未知量s的形式：
|
(11) |
其中，d是已知的原始地震信号，s为未知的有效信号成分，如何准确获取有效信号是局部信噪比计算的关键.
2.2 有效信号估计Chen和Fomel (2015)利用信号和噪声正交的原理提出了一种正交更新方法提取信噪分离后噪声中残存的信号.初始信噪分离后的分离出信号s0和噪声中残留的信号s1间满足如下线性关系：
|
(12) |
其中，ω是正交系数，易知diag (ω)s0=diag (s0)ω.
真实有效信号 (或噪声) 和信噪分离后信号 (或噪声) 的关系为
|
(13) |
|
(14) |
其中，s是恢复后的真实有效信号，n是去除了有效信号的噪声 (或称之为真实噪声)，n0是初始信噪分离出来的噪声，s0+n0=d.假设s和n是正交的，则
|
(15) |
把式 (13) 和 (14) 代入式 (15)，可以计算出正交系数ω为
|
(16) |
根据上述推导可知，求取正交系数的前提条件是有效信号s和真实噪声n正交.然而当出现初始信噪分离后信号中残留噪声的情况 (简称为“欠滤波”) 时，该方法是在求含噪信号与真实噪声的正交解，违反了前提条件.此时式 (13) 和式 (14) 不再适用，可以得到如下新的正交关系式：
|
(17) |
其中，n1是信噪分离后信号中残留的噪声.但是n1的引入导致公式 (17) 成为欠定问题，造成方程无法求解.为了解决这个问题，本次研究提出基于“过滤波”(初始滤波后噪声剖面中包含有效信号，而滤波后信号中无噪声) 的级联信号估计方法.当保证“过滤波”有效时，可以获得准确的正交算子.而相对“欠滤波”而言，在应用很多常规滤波器时，如平滑滤波器，可以达到预期的“过滤波”效果.该方法相当于两个滤波器的级联组合，通过第一阶段的二维过滤波器使滤波后的信号中几乎不含噪声，以此满足第二阶段正交更新的前提条件，从而求得正确的正交系数.
多维高斯滤波器是一种常用的滤波方法，其具有较强的去噪能力和较差的信号保护能力.实际中该滤波器除了能去除高斯分布的随机噪声，对其他类型的随机噪声也具有较强的去噪能力，因此本次研究选取其作为过滤波器.
对三角窗函数进行傅里叶变换得到窗谱函数为
|
(18) |
其中，f是信号频率，T是窗长度.公式 (18) 的p阶卷积为
|
(19) |
根据中心极限定理可知，当p足够大时，公式 (19) 所示的多阶褶积滤波器收敛为高斯滤波器，即通过对三角平滑滤波器的自身褶积可以实现快速的高斯滤波器.
高斯过滤波后的信号为
|
(20) |
其中，w(t, x) 是过滤波器的时空域算子.这样，过滤波后的信号就可以代入公式 (13) 和公式 (14) 来表示正交系数.得到正交系数的表达式后，不难发现公式 (16) 和公式 (5) 具有相似形式，两者都可转化为对不适定反问题的求解，因此公式 (16) 可以写为如下正则化最小二乘形式：
|
(21) |
其中，R是正则化平滑算子.Fomel (2005)提出一种整形正则化方法，该方法具有较小的条件数和较快的收敛速度，而且具有自适应变化的平滑窗口.整形正则化条件约束下正交系数ω的最小二乘解为
|
(22) |
其中，S0=diag (s0)，Γ是整形算子，λ是尺度系数，λ=‖S0TS0‖2.把求出的正交系数代入式 (13) 得到真实有效信号为
|
(23) |
将公式 (23) 代入公式 (11)，得到最终的局部信噪比计算公式.
3 模型测试在理论模型测试中，使用一个包含倾斜平行同相轴、正弦形态同相轴、断层和不整合面的叠后模型数据 (Claerbout, 2009) 进行测试，如图 1a所示.加入随机噪声的数据如图 1b.将理论有效信号代入公式 (2) 求出含噪信号的全局信噪比为3.128，再把有效信号分别代入公式 (3) 和公式 (11) 可以求出含噪模型的瞬时信噪比 (图 1c) 和局部信噪比 (图 2).通过图 1c可以看出瞬时信噪比剖面存在很多异常点，无法正确反映信号的信噪比关系.图 2所示的局部信噪比则很好地对应了全局信噪比数值，而且刻画了信号的局部信噪比特征，证明了局部信噪比的表示方法是正确的.图 2给出了三种不同平滑半径为参数的正则化约束下信噪比图.平滑半径越小，得到的信噪比越精确，随着平滑半径增大，局部信噪比的分辨率降低，但是受零点异常的影响减小，方法越稳健，因此在实际应用中应视具体情况来选择平滑半径.接下来对含噪信号使用级联信号估计方法计算有效信号.选取三角滤波器平滑半径为3，迭代次数为10，将公式 (23) 的计算结果代入公式 (2) 得到全局信噪比为3.672 dB，进一步得到含噪模型的局部信噪比如图 3所示.对比图 2和图 3可以看出两者已比较接近.全局信噪比和局部信噪比的双重对比证明级联信号估计方法能够在有效去除随机噪声的同时最大程度地避免有效信号的损失，使求出的局部信噪比更加接近真实值.最后，使用工业标准的f-x预测滤波 (Liu et al., 2015) 表明局部信噪比对评价去噪方法的作用.f-x预测滤波具有对近似线性同相轴滤波效果较好，对正弦形态同相轴滤波效果较差的特点.对含噪信号模型进行f-x预测滤波，去噪结果如图 4a所示，求取去噪后信号的全局信噪比为9.972 dB.接下来对去噪后的有效信号进行局部信噪估计，设置平滑半径为10，结果如图 4b所示.可以看出图 4b的数值范围和全局信噪比数值有较好的对应.对滤波前后的局部信噪比作差并进行归一化，结果如图 4c所示.可以看出滤波后地震信号信噪比在倾斜平行同相轴处得到明显提高，在正弦形态同相轴处信噪比变化较小，这和f-x预测滤波的滤波特点得到了较好的对应.另外，在倾斜平行同相轴的上层部分地震信号的信噪比提升不大，可能为f-x预测滤波在倾斜平行同相轴处产生虚假的正弦形态同相轴所产生的影响.理论模型测试表明，局部信噪比方法不仅能够直观准确地指示地震信号质量好坏和局部信噪比特征，而且对去噪方法的定量评价有着重要意义.
|
图 1 含随机噪声模型及其瞬时信噪比 (a) 真实有效信号；(b) 含噪信号；(c) 瞬时信噪比. Fig. 1 Noisy model and its instantaneous SNR (a) Ideal signal; (b) Noisy signal; (c) Instantaneous SNR. |
|
图 2 利用真实有效信号求出的含噪信号局部信噪比 (a) 平滑半径为5；(b) 平滑半径为10；(c) 平滑半径为15. Fig. 2 Local SNR for noisy data by using ideal signal (a) Smooth radius is 5; (b) Smooth radius is 10; (c) Smooth radius is 15. |
|
图 3 级联信号估值方法求出的含噪信号局部信噪比 (a) 平滑半径为5；(b) 平滑半径为10；(c) 平滑半径为15. Fig. 3 Local SNR for noisy data by using cascading signal estimation (a) Smooth radius is 5; (b) Smooth radius is 10; (c) Smooth radius is 15. |
|
图 4 f-x预测滤波后的地震信号及其局部信噪比 (a) f-x预测滤波后的地震信号；(b) 滤波后局部信噪比 (平滑半径为10)；(c) 归一化的滤波前后局部信噪比差值. Fig. 4 Data processed by using f-x deconvolution and it local SNR (a) Signal estimated by f-x deconvolution; (b) Local SNR (smooth radius is 10); (c) Local SNR difference between before and after filtering. |
为了进一步测试本次研究方法的有效性，使用如图 5所示的二维叠后实际地震资料 (Liu and Chen, 2013) 进行局部信噪比估计.该地震资料共有310道地震记录，时间采样间隔为1 ms，其中的随机噪声主要是由该地区的地表条件所引起.首先，对原始地震剖面使用级联信号估计方法计算信号成分，分别给出过滤波后的有效信号 (图 6a) 和噪声 (图 7a) 及正交恢复后的有效信号 (图 6b) 和噪声 (图 7b).为证明级联信号估计方法中过滤波器的必要性，给出欠滤波方法与正交更新方法级联的结果.选取0.1~50 Hz的频率域带通滤波器作为欠滤波器，滤波结果如图 6c所示.正交更新方法级联后计算出噪声剖面如图 7c所示.对信号和噪声做相关性分析 (Fomel, 2007)，如图 8所示，相关性数值越大，说明噪声和有效信号的相关性越强，也就是噪声中含有的有效信号越多.可以看出，级联信号估计方法计算的相关性远低于过滤波方法和欠滤波正交更新方法.将估计出的有效信号代入信噪比公式分别计算原始地震信号的全局信噪比 (15.106 dB)、瞬时信噪比 (图 9a) 和局部信噪比 (图 9b).图 9a所示的瞬时信噪比数值范围远超全局信噪比，没有参考意义.图 9b是平滑半径为10的局部信噪比，其数值范围和全局信噪比有较好对应.在图 9b中，方框处是信号局部信噪比较高的地方，而在图 5的相同位置，原始地震剖面的有效波同相轴振幅较强，两者有较好的对应关系.由此可见，局部信噪比估值方法能够有效反映地震资料的局部信噪比特征，从而可靠地判断地震资料的质量好坏.另外为了说明局部信噪比估计方法对去噪方法的评估作用，对原始地震剖面使用f-x域预测滤波去噪，去噪后的地震信号如图 10a所示.使用局部信噪比估计方法对图 10a进行计算得到如图 10b所示的局部信噪比图像，对比图 9b和图 10b，可以看出f-x预测滤波去噪后地震信号的信噪比得到了提高，尤其是在线性同相轴区域去噪效果明显，但在弯曲同相轴处去噪能力较差.对滤波前后的局部信噪比作差并进行归一化，结果如图 10c.可以看到线性同相轴处的信噪比都得到了明显提高，弯曲同相轴处的信噪比提高程度则相对较低.0.1~0.25 s之间的线性同相轴在滤波前信噪比为25 dB左右，远高于其他部分的信噪比，使得归一化时分母基数较大，所以在图 10c中该处数值较低.
|
图 5 二维地震资料 Fig. 5 Two-dimension field data |
|
图 6 处理后的地震资料 (a) 过滤波后地震信号剖面；(b) 级联估计后地震信号剖面；(c) 欠滤波后地震信号剖面. Fig. 6 Signal after different processing (a) Over filtering; (b) Cascading estimation; (c) Insufficient filtering. |
|
图 7 地震资料噪声剖面 (a) 过滤波后噪声剖面；(b) 级联估计后噪声剖面；(c) 欠滤波正交更新后噪声剖面. Fig. 7 Removed noise after different processing (a) Over filtering; (b) Cascading estimation; (c) Insufficient filtering. |
|
图 8 信号和噪声相关性 (a) 过滤波后信号和噪声相关性；(b) 级联估计后信号和噪声相关性；(c) 欠滤波正交更新后信号和噪声相关性. Fig. 8 Correlation between signal and noise after different processing (a) Over filtering; (b) Cascading estimation; (c) Insufficient filtering. |
|
图 9 地震资料局部信噪比 (a) 瞬时信噪比；(b) 局部信噪比 (平滑半径为10). Fig. 9 SNR for the field data (a) Instantaneous SNR; (b) Local SNR (smooth radius is 10). |
|
图 10 f-x滤波后地震资料及其局部信噪比 (a) f-x滤波后地震信号剖面；(b) f-x滤波后地震信号剖面的局部信噪比；(c) 归一化的滤波前后局部信噪比差值. Fig. 10 Data after f-x deconvolution and its local SNR (a) Data filtered by using f-x deconvolution; (b) Local SNR for filtered data; (c) Local SNR difference before and after f-x deconvolution. |
实际资料试算结果表明，局部信噪比估计方法可以准确完成信噪比的定量分析和去噪方法的定量评价，为进一步的处理工作提供可靠依据.
5 结论非平稳地震信号的自适应处理有着十分重要的意义，然而目前为止鲜有针对非平稳信号信噪比的局部估计方法.鉴于此，本文提出了一种基于正则化条件的局部信噪比方法.针对全局信噪比无法表征信号的局部特征而瞬时信噪比计算结果不准确的问题，局部信噪比估计方法利用信号中每个数据点及其邻域各点的局部信息，避免因为使用单个数据点信息而可能出现的不合理值.此外，本文改进了一种级联信号估计方法以避免已有的正交更新方法的失败情况.理论模型测试和实际资料处理结果均表明，局部信噪比估值方法能够准确反映地震资料的局部信噪比特征，为判断地震资料质量的好坏以及评估去噪方法的定量关系提供了可靠依据.
致谢感谢中石油东方地球物理公司物探技术研发中心罗国安专家提供的技术思路和讨论.
| Álvarez-DíazM, LóPez-ValcarceR, MosqueraC. 2010. SNR estimation for multilevel constellations using higher-order moments. IEEE Transactions on Signal Processing, 58(3): 1515–1526. | |
| Benedict T R, Soong T T. 1967. The joint estimation of signal and noise from the sum envelope. IEEE Transactions on Information Theory, 13(3): 447-454. DOI:10.1109/TIT.1967.1054037 | |
| Chen Y K, Fomel S. 2015. Random noise attenuation using local signal-and-noise orthogonalization. Geophysics, 80(6): WD1-WD9. DOI:10.1190/geo2014-0227.1 | |
| Claerbout J F. Basic Earth Imaging: Stanford Exploration Project. California: Stanford University Press, 2009. | |
| Fomel S. 2007. Shaping regularization in geophysical-estimation problems. Geophysics, 72(2): R29-R36. DOI:10.1190/1.2433716 | |
| Fomel S. 2007. Local seismic attributes. Geophysics, 72(3): A29-A33. DOI:10.1190/1.2437573 | |
| Huo Z J, Fu D D. 1990. Estimating quantitatively S/N of seismic record. Journal of Jianghan Petroleum Institute, 12(3): 15-21. | |
| Kerr R B. 1966. On signal and noise level estimation in a coherent PCM channel. IEEE Transactions on Aerospace & Electronic Systems, AES-2(4): 450-454. | |
| 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. DOI:10.1190/geo2014-0011.1 | |
| Liu G C, Chen X H. 2013. Noncausal f-x-y regularized nonstationary prediction filtering for random noise attenuation on 3D seismic data. Journal of Applied Geophysics, 93: 60-66. DOI:10.1016/j.jappgeo.2013.03.007 | |
| Liu Y, Li C C. 1997. Some methods for estimating the signal/noise ratio of seismic data. OGP, 32(2): 257-262. | |
| López-Valcarce R, Mosquera C. 2007. Sixth-order statistics-based non-data-aided SNR estimation. IEEE Communications Letters, 11(4): 351-353. DOI:10.1109/LCOM.2007.348298 | |
| Matzner R. 1993. An SNR estimation algorithm for complex baseband signals using higher order statistics. Facta Universitatis (Nis), 6(1): 41-52. | |
| Niu C, Zhan Y, Li H F. 2006. The contrast of several methods for the signal/noise ratio estimation of seismic records. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 28(1): 5-9. | |
| Pratt W K. Digital Image Processing. New York: Wiley, 1978. | |
| Tikhonov A N, Arsenin V Y. Solutions of Ill-Posed Problems. New York: Winston, 1977. | |
| Turin G. 1960. An introduction to matched filters. IRE Transactions on Information Theory, 6(3): 311-329. DOI:10.1109/TIT.1960.1057571 | |
| Wang H L. 2007. The research of methods to estimate signal-to-noise ratio of seismic record[Ph.D.thesis] (in Chinese). Chengdu: Chengdu University of Technology. | |
| Zhang J H, Zang S T, Zhou Z X, et al. 2009. Quantitative computation and comparison of S/N ratio in seismic data. OGP (in Chinese), 44(4): 481-486. | |
| Zhao L, Hu Z F, Song J P, et al. 2003. Comprehensive evaluation of seismic data attribute. Henan Petroleum (in Chinese), 17(1): 6-8. | |
| 霍志坚, 傅旦丹. 1990. 地震记录信噪比的定量估算. 江汉石油学院学报, 12(3): 15–21. | |
| 刘洋, 李承楚. 1997. 地震资料信噪比估计的几种方法. 石油地球物理勘探, 32(2): 257–262. | |
| 牛聪, 詹毅, 李辉峰. 2006. 对比地震记录信噪比的几种估算方法. 物探化探计算技术, 28(1): 5–9. | |
| 王红玲. 2007. 地震记录信噪比估算方法研究[博士论文]. 成都: 成都理工大学. | |
| 张军华, 藏胜涛, 周振晓, 等. 2009. 地震资料信噪比定量计算及方法比较. 石油地球物理勘探, 44(4): 481–486. | |
| 赵亮, 胡志方, 宋建平, 等. 2003. 地震资料品质综合评价方法. 河南石油, 17(1): 6–8. | |