随机噪声的局部正交压制方法

徐彦凯 曹思远* 潘 晓 杨光明 张宪旭

(①中国石油大学(北京)地球物理与信息工程学院，北京 102249； ②中国石油大学(北京)CNPC物探重点实验室，北京 102249；③南方科技大学地球与空间科学系，广东深圳 518055； ④中煤科工集团西安研究院有限公司，陕西西安 710077)

0 引言

噪声一直是影响地震资料品质的重要因素之一，地震资料处理的最基本目标就是消除地震记录中的各种噪声，最大程度地提高地震资料信噪比。对随机噪声的压制许多学者已提出了大量的方法并取得了很多研究成果，这些方法主要有中值滤波方法[1]、f-x反褶积方法[2]、离散余弦变换方法[3]、小波变换方法[4-6]、独立分量分析方法[7]、经验模态分解[8]、f-x预测滤波方法[9-10]、字典学习方法[11]以及奇异值分解(SVD)[12-15]方法等。

刘财等[1]将二维多级中值滤波应用到地震勘探中的噪声压制，研究了噪声强度与滤波因子的关系。陆文凯[3]在离散余弦变换域中利用预测滤波器压制噪声。Cao等[5]给出了一种柔性的基于提升法的小波构造方法，适用线性、非线性或空间变化的预测和更新算子，并能确保变换的可逆性。刘洋等[6]采用高阶seislet变换的提升算法对地震信号去噪，重点考虑了去噪阈值选取。Bekara等[8]在f-x域采用经验模态分解降噪； 以此为基础，Chen等[9]在f-x域中先使用倾角滤波器，再利用经验模态分解压制了复杂地层结构中的噪声。Liu等[10]在f-x域中采用整形正则化约束预测倾角压制噪声，滤波器长度和整形算子是该方法中二个重要参数。陆文凯等[12-13]采用SVD方法实现地震资料的降噪，之后在纹理方向检测倾角的基础上再用SVD方法压制噪声。徐彦凯等[15]采用奇异值选择的自适应方法压制地震资料中的噪声。

上述方法均探讨了如何将信号与噪声有效分离。当然，信噪分离效果的好坏与压制噪声的阈值、滤波器算子长度的大小有很大关系。阈值较小或滤波器算子较短，信噪分离的程度相对较小，即去噪后的地震资料中还存在较多噪声； 反之，若将噪声完全从地震资料中去除，地震信号可能会受到部分损伤，即去除的噪声中尚存部分有效信号。为了提高地震资料信噪比，本文充分利用了信号和随机噪声的不相关性，依据信号与噪声满足的局部正交化准则，从分离后的噪声中再提取地震信息。

1 基本原理

地震数据d可看成由理想地震信号s和随机噪声n两部分组成，即

d=s+n≈s0+n0

(1)

式中s0和n0分别为采用传统方法处理得到的信号与噪声，即第一次分离出的信号与噪声。

大多数情况下，第一次分离噪声n0不等于n和n0中含有的残留信号s1。从图1a～图1d、图1f、图1g可见：两种方法分离的噪声中都含有一些信号的影子。对分离后的信号与噪声做相关分析，得到相关系数剖面(图1e、图1h)，可见分离出的噪声与信号有很强相关性。依据其相关特性估计噪声中的残留信号。残留信号s1由分离出的信号s0加权[16]表示为

s1=w○s0=Ws0=S0w

(2)

式中：w表示通过s0估计s1的权重因子； “○”表示hadamard积；W=diag(w)；S0=diag(s0)。

1.1 局部正交准则

=s0+ws0

(3)

(4)

式(3)和式(4)中w表示全局权重系数。

由于信号s与随机噪声n不相关，则与相互正交，即

(5)

将式(3)和式(4)代入式(5)，可得

(6)

图1 常规处理后信号与噪声及相关性

由此得到全局权重系数

(7)

加窗后式(7)改写为[16]

(8)

w(t)是t点处窗长m区域内满足局部正交准则的权重系数。式(8)的矩阵形式为

(9)

与全局正交相比，采用局部正交准则能根据有效信号缺失程度自适应地优化权重因子w。

从式(9)可见：当已知S0和n0时，计算权重因子w是不适定问题。对此，可采用整形正则化求解。

1.2 整形正则化

正则化目的是对估计的模型解强加限制，以使不适定问题得以求解，其中Tikhonov正则化为一种经典方法。在此基础上，Fomel[17]考虑整形算子的平滑作用，提出了整形正则化理论。

式(9)可转化为最小平方意义下目标函数优化

(10)

对目标函数加约束，式(10)改写为

(11)

式中：D为对权重因子w加以光滑约束的Tikhonov 正则化项；ε为正则化系数。采用最小二乘法得到其解为

(12)

求解式(11)的优化问题时，考虑解的平滑性和收敛速度，引入约束模型在可接受空间映射的整形算子[18]

Γ=(I+ε2DTD)-1

(13)

式(12)可改写为

(14)

将式(14)代入式(3)，得到信号s的估计。

1.3 自适应奇异值分解(ASVD)

设地震道记录中每道有N个采样点。该记录可用矩阵

(15)

表示，其元素为aij(i为道号，j为时间样点号)。

设矩阵AN×M(N≤M)的秩为r，则存在两正交矩阵U、V和对角矩阵Λ，且满足

AN×M=UΛVT

(16)

其中

式中：Δ=diag[σ1,σ2,…,σr],σ1≥σ2…≥σr>0；σi是矩阵A的奇异值；ui、vi分别是矩阵A的左、右奇异向量。

地震信号经过奇异值分解后，奇异值仅集中在前几个。这样，信号奇异值分布有一个陡降，而随机噪声奇异值由大到小的递减幅度统计近似相同，因此采用相邻奇异值增量分离地震资料中的随机噪声。相邻奇异值增量定义为

Δσi=σi-σi+1

(17)

ASVD是将数据做SVD得到奇异值后，计算奇异值增量的幅度和序号以及阈值，在序号内确定最大奇异值增量，由阈值确定保留奇异值的个数，从而实现信号和噪声分离。

2 模型与实际资料处理

2.1 模型试验

为了检验本文算法的有效性，首先对图1b做ASVD，图2a与2b分别为分离得到的信号与噪声，可见图2a中噪声很少； 计算二者的相关性(图2c)后与图1e和图1h比较，可见由ASVD得到的信号与噪声的相关性更小。

在此基础上，采用整形正则化的局部正交法从分离出的噪声中提取残留信号并进行处理，然后与分离前的相关系数(图1e、图1h)及采用f-x[16]、常规SVD方法的结果进行对比(图2)。可见分离后信号(图2d、图2g、图2j)与噪声(图2e、图2h、图2k)的相关系数(图2f、图2i、图2l)都有不同程度的减小； 本文方法得到的有效信号剖面更干净清晰，相关系数明显小于文献方法的结果[16]。

图3a、图3c、图3e分别为理想信号与用f-x、常规SVD和ASVD三种方法压制噪声后信号的误差剖面； 图3b、图3d、图3f为分别为f-x、常规SVD和ASVD三种方法基础上采用局部正交法分离得到的误差剖面。对比这两组剖面可见：f-x和ASVD两种方法采用局部正交法处理后的误差小于未采用的；采用局部正交法分离将f-x方法的误差剖面(图3a)上水平同相轴信息得到了充分提取(图3b)；由于ASVD方法中倾斜同相轴不能全部集中在前几个奇异值上，所以图3e中倾斜同相轴误差较大，采用局部正交法分离将ASVD方法的误差剖面(图3e)中倾斜同相轴信息实现了较好地提取，得到误差剖面图3f。

抽取图1中第10道的理想信号(图4a)和加噪信号(图4b)，以及f-x、常规SVD和ASVD三种方法处理结果(图4c、图4e、图4g)、对应的局部正交法处理后结果的第10道(图4d、图4f、图4h)，可见用ASVD方法+局部正交法(图4h)压制噪声效果最好。为了便于量化对比几种方法的压噪效果，计算了局部正交法分离前、后的信噪比(表1)。

表1 局部正交法分离前、后信噪比

图2 局部正交法压噪效果

图3 图2对应的误差剖面

图4 图2中第10道压噪效果

表1中f-x和ASVD分离信噪比较高，信号中含随机噪声较少，通过局部正交法分离的噪声中能有效提取残留信号，从而提高分离后的信噪比。 SVD方法分离的信号中噪声较多，采用局部正交法从分离的噪声提取的残留信号中噪声较多，因此分离后的信噪比反而减小。由上可知，采用局部正交分离要求分离前的信号含噪少。

为了进一步验证局部正交法信噪分离效果，对水平地层模型做ASVD降噪处理[15]后，又采用局部正交法提取分离噪声中的残留信息。设炮点到各层的最小时间t01～t05分别为0.1、0.4、0.6、0.8、0.95s，对应速度v1～v5分别为950、1000、1100、1200、1500m/s； 采样间隔为4ms； 采样点数为300；道间距为10m； 道数为81； 雷克子波峰值频率为35Hz。在合成的CDP道集(图5a)上，加入一定噪声后得到信噪比为-1dB的记录(图5b)，并对其进行无拉伸动校处理(图5c)； 应用ASVD进行信噪分离，得到一次分离后的信号(图5d)、噪声(图5e)及其二者的相关系数(图5f)； 采用局部正交法从一次分离的噪声中提取残留信号，得到最终的信号(图5g)、噪声(图5h)及其二者的相关系数(图5i)； 对图5d和图5g分别做反无拉伸动校处理，得到对应地震记录(图5j、图5k，其信噪比分别为6.4、7.5dB)、一次分离后(图5l)及最终(图5m)噪声。从二次分离前后的信号、噪声及其相关性可见：局部正交法可有效提取噪声中的残留信号，大大降低信号与噪声的相关性，从而最大限度地保留有效信号并压制噪声，进而提高地震记录信噪比。

(a)合成CDP道集； (b)信噪比为-1dB的含噪记录； (c)无拉伸动校后记录； (d)、(e)、(f)分别为ASVD第一次分离后信号、噪声及二者相关系数； (g)、(h)、(i)局部正交法处理后的信号、噪声及二者相关系数； (j)、(k)分别是图d、图g的反无拉伸动校后记录； (l)一次分离后噪声；(m)最终噪声

2.2 实际资料处理

选取实际二维叠加数据(图6)的第300～500道，其采样间隔为2ms，每道400个采样点。进行f-x处理(图7a)和局部正交处理(图7c)，得到相应的差剖面(图7b、图7d)；再做ASVD处理(图7e)和局部正交处理(图7g)，得到对应的差剖面(图7f、图7h)。从上述差剖面中方框部分可见： 局部正交处理前的差剖面(图7b、图7f)都含有明显的同相轴信息，而局部正交处理后的差剖面(图7d、图7h)中同相轴信息较好地被提取；对比局部正交处理后的f-x方法(图7d)与ASVD方法(图7h)差剖面，可见后者处理效果更好； 同时，其对应的叠加剖面图7g的信噪比高于图7c。

图7 实际地震数据的不同处理方法剖面对比(a)、(c)分别为f-x和f-x+局部正交处理的叠加剖面； (b)、(d)分别为f-x和f-x+局部正交处理的差剖面； (e)、(g)

3 结束语

地震数据处理中，常规方法压制随机噪声后，滤除噪声中尚存一定量有效信号。本文提出的基于整形正则化的局部正交方法可从常规分离的噪声中有效提取残留信号。该方法要求常规分离的信号中噪声较少。与f-x、SVD方法相比，ASVD方法压噪性能更强。模型和实际资料处理结果均表明，在ASVD降噪基础上采用局部正交法可有效压制地震数据中随机噪声，显著提高地震资料信噪比。