地球物理学报  2021, Vol. 64 Issue (1): 263-278   PDF    
基于Contourlet域自适应Wiener阈值的同时震源波场分离
王坤喜1,2, 毛伟建1     
1. 中国科学院精密测量科学与技术创新研究院计算与勘探地球物理研究中心; 大地测量与地球动力学国家重点实验室, 武汉 430077;
2. 中国科学院大学, 北京 100049
摘要:同时震源数据包含了多炮之间的串扰噪声,不能直接用于常规数据处理流程.因此,需要对混叠的波场进行分离得到常规采集的单炮记录.本文基于稀疏迭代反演分离,提出了一种具有尺度与空间自适应的Wiener阈值选取方法.该阈值选取方法能够根据不同迭代环境计算不同尺度下串扰噪声的方差和不同空间位置有效信号的方差,从而自适应调整阈值大小,最终通过对变换域系数进行收缩来达到去除串扰噪声的目的.理论模型数据和实际数据测试结果表明,本文方法能够快速有效地压制串扰噪声和保护弱有效信号,取得了比Contourlet域子带一致Wiener阈值方法和Curvelet域指数衰减阈值方法更好的分离效果.
关键词: 同时震源      稀疏反演      Contourlet变换      自适应阈值      Wiener滤波     
Simultaneous-source wavefield separation based on adaptive Wiener threshold in Contourlet domain
WANG KunXi1,2, MAO WeiJian1     
1. Center for Computational and Exploration Geophysics, Innovation Academy for Precision Measurement Science and Technology, Chinese Academy of Sciences; and State Key Laboratory of Geodesy and Earth's Dynamics, Wuhan 430077, China;
2. University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: The simultaneous-source data cannot be directly used for seismic data processing workflow due to crosstalk noise of the interference guns. Therefore, the blended wavefield needs to be separated to obtain the conventional single shot record. Based on the threshold iterative inversion separation in sparse domain, a scale and spatially adaptive Wiener threshold method is proposed. With the method, the threshold values are adaptively calculated by using the variance of crosstalk noise at different scales and the variance of effective signals at different spatial positions at each iteration, and are used to attenuate the crosstalk noise by shrinking the coefficients in transform domain. The results tested by synthetic and field data show that the proposed method can more quickly and effectively suppress crosstalk noises and keep valuable weak signals than the traditional Wiener threshold method with subband consistent in Contourlet domain and the exponential attenuation threshold method in Curvelet domain.
Keywords: Simultaneous-source    Sparse inversion    Contourlet transform    Adaptive threshold    Wiener filtering    
0 引言

同时震源采集技术,也称多震源地震数据混合采集技术,可以在常规采集一个单炮地震记录的时间内采集多个炮的地震记录.当施工的周期不变时,通过增加激发的震源数量,提高对地下构造的覆盖次数,改善对地下构造的照明质量;当对目标区域的覆盖次数不变时,多个分布在不同区域的震源同时激发,可以大大缩短施工周期,提高采集效率(Li et al., 2013Zu et al., 2017张攀和毛伟建,2018).由于相邻震源之间的干扰会产生“混合噪声”,同时震源数据不能直接用于传统地震数据处理流程,需要先将混叠波场进行分离处理,得到常规的单炮数据.分离方法通常可以分为被动分离和主动分离两大类(Akerberg et al., 2008Liu et al., 2014Magnussen,2015).

被动分离是将分离问题转换为去噪问题,将经过延迟时间编码的干扰炮信号视为噪声.当不同震源之间存在随机时间差时,混合噪声在非共炮域(如共检波点域、共中心点域、共偏移距域等)不再具有连续性,表现为伪随机脉冲,而有效信号仍然保持连续相干性(Doulgeris et al., 2012).基于上述性质,Moore等(2008)在共偏移距道集进行滤波分离.Huo等(2009)利用多方向矢量中值滤波法对地层的倾角进行了扫描,求得最佳倾角方向后利用中值滤波去除共中心点域的混合噪声.Kim等(2009)基于数据本身建立了一个噪声估计模型,然后从采集数据中自适应的提取噪声,该算法适用于共炮检距道集.Zhang和Olofsson(2012)采用加权τ-p变换压制邻炮干扰.Yang等(2017)提出先用PWD(plane-wave destruction)求取地层的倾角信息,再结合矢量中值滤波对串扰噪声进行压制.

另一种分离类方法是主动分离,也被称为反演类方法,该方法是将分离问题提炼成数学中的反问题,利用反演类方法进行求解.Akerberg等(2008)Moore(2010)刘强等(2014)在Radon域进行了稀疏约束下的迭代反演分离.Lin和Herrmann(2009)在Curvelet域进行稀疏表示,利用压缩感知对混合的波场进行分离.Abma等(2010)在频率域通过稀疏约束的POCS(projection onto convex sets)方法使不相干信号能量最小化,进而获取分离信号.Tan等(2013)基于互易定理,利用稀疏反演的方法进行波场分离.Chen等(2014)提出整形正则化反演的方法并在Seislet域利用阈值进行约束反演.祖绍环等(2016)在Curvelet域进行稀疏约束,同样利用阈值进行约束反演.宋家文等(2019)将地震数据转换到频率—波数—波数域(FKK),提出了一种高效混采数据分离的方法.在稀疏迭代分离框架下,除了使用固定的变换域外,一些学者也考虑采用字典学习来完成分离(Chen,2017Zhou et al., 2017),但是字典学习的方法需要很大的计算开销.

在同时震源稀疏域迭代反演中,选择一个合适的稀疏基十分重要, 它决定了分离效果和迭代收敛的速度.Do和Vetterli(2002)首次提出了Contourlet变换理论,将拉普拉斯金字塔(Burt and Adelson, 1983Do and Vetterli, 2003)和二维方向滤波器组串联在一起组成金字塔多方向滤波器组(Do and Vetterli, 2001),可以将信号数据分解成不同尺度下的方向子带,很好地抓住数据的几何结构特征,满足曲线的各向异性尺度关系,能够快速并结构化分解采样信号.Contourlet变换将多尺度分析和方向分析分开进行,首先由拉普拉斯金字塔变换对图像进行多尺度分解以“捕获”点奇异,接着由方向滤波器组将分布在同方向上的奇异点合成一个系数.Contourlet变换的最终结果是用类似于线段的基结构来逼近原图像,相对于Wavelet变换和Curvelet变换等多尺度分析工具,Contourlet变换具有更少的冗余度和更好的逼近性能,因此Contourlet变换已在信号处理、图像处理等领域得到广泛应用,如信号去噪(Eslami and Radha, 2003),目标检测(Wilbur et al., 2009),图像压缩(Chang et al., 2000a)等等.

与数字图像去噪的目标相似,同时震源的稀疏域迭代反演分离需要在去除串扰噪声的同时尽可能地保留地震数据中的弱信号和细节信息(Zu et al., 2017),合理地选取阈值非常重要,直接影响到分离效果(Yang et al., 2013).Donoho和Johnstone(1994)在小波域提出了阈值收缩方法,给出了Universal阈值,并从渐进意义上证明了Universal阈值的最优性(Donoho,1995),然而Universal阈值是阈值的上限,但不是最佳收缩阈值,它过度“扼杀”了小波系数,丢失信号的细节部分.Donoho和Johnstone(1995)又提出了SURE阈值,但是该阈值偏于保守,在信噪比很小的情况下,其效果不如Universal阈值.Zhang和Desai(1998)改进SURE阈值提出了SureShrink阈值,可以求出理想阈值的估计值,但没有显式表达式,阈值的计算往往需要预先知道信号本身,但实际求取中这是不可能的.Chang等(2000a)假设原始不带噪声信号在小波域的系数服从广义高斯分布,在贝叶斯框架下通过最小化贝叶斯风险,得到了著名的BayesShrink阈值,该阈值考虑了信号在小波域的先验知识,具有很好的去噪效果.在借鉴图像压缩编码中的上下文模型方法后,Chang等(2000b)又提出了具有局部自适应能力的BayesShrink阈值,能够灵活实现去噪与保留信号之间的平衡.在贝叶斯最大后验概率估计理论下,Moulin和Liu(1999)假设图像子带的小波系数服从拉普拉斯分布或者高斯分布,可以分别推导得到MapShrink阈值与Wiener阈值,这两个阈值算子都具有较好的去噪效果.

以上阈值方法的提出主要用于数字图像去噪中,一些学者也在研究适用于同时震源稀疏域迭代反演分离的阈值算子.Chen等(2014)提出在Seislet域采用百分位法(Yang et al., 2013)分段设置硬阈值,取得了一定的分离效果,后来又在PNMO-MF-FK联合滤波器下使用线性衰减硬阈值也得到了较好的分离结果.该阈值算子在初始迭代时取系数最大值,并随着迭代次数线性衰减为零(Chen et al., 2015).Li等(2018)将混叠数据变换到Curvelet域后,发现相干信号主要集中在低频子带中,串扰噪声主要集中在高频子带中,通过考虑在不同子带设置不同的阈值,取得了比全局一致方法更好的效果.Zu等(2018)将该阈值方法推广并运用于三维同时震源数据分离.在同时震源稀疏域迭代反演分离方法中,另一种常见的阈值是指数衰减阈值(Gao et al., 2010Cheng and Sacchi, 2015Song et al., 2019).与常数和线性衰减阈值相比,指数衰减阈值具有更快的收敛速度,但是该阈值没有充分利用稀疏域系数的分布情况等先验信息,不具备自适应去噪的能力.

本文利用Contourlet变换良好的稀疏表达能力,通过稀疏域阈值收缩来迭代反演分离同时震源数据.我们在子带一致的Wiener阈值基础上做出改进,提出具有尺度与空间自适应的Wiener阈值.通过两个模拟数据和两个实际数据的应用结果验证本文提出的自适应Wiener阈值方法的有效性.

1 方法 1.1 稀疏迭代反演分离

在时间域,每个检波点接收到的同时震源混合记录可表示为

(1)

其中dobs为混叠道集记录,db为第b个震源的原始道集记录,m为混叠度.对角矩阵Γb表示第b个震源的延迟时间矩阵,其表达式为

(2)

其中,tb, u为第b个震源第u道激发的延迟时间,i为虚数单位,ω为角频率,表示一对正逆傅里叶变换.通常,在同时震源数据分离中,第一个震源(本文称为主炮)的延迟时间算子(Γ1)为单位矩阵I,得到连续相干的有效信号;其他震源(本文统称为干扰炮)在延迟时间算子的作用下会形成串扰噪声.

同时震源数据的稀疏反演分离问题可以表示如下:

(3)

其中‖·‖22表示2范数.式(3)是欠定的(Moore,2010韩立国等,2013魏亚杰等,2019),需要额外的约束来获得唯一的解.

求解式(3)的一个常用迭代框架是非线性整形正则化框架(Chen et al., 2014):

(4)

其中,n为迭代次数,Γg-1表示第g个震源对应的延迟时间算子的逆.S为整形算子,对模型具有约束作用,其表达式为

(5)

其中T为阈值参数,CC-1分别代表稀疏变换及其逆变换,本文中C为Contourlet变换.

共检波点道集或共偏移距道集中的主炮信号在Contourlet域是聚焦的能量点,而干扰炮信号是分散的串扰噪声,T可以对变换域系数进行收缩从而有效地压制串扰噪声进而分离混叠波场.因此,同时震源数据的稀疏迭代反演分离的关键是如何确定合适的T.

1.2 子带一致Wiener阈值

在每一次迭代下使用阈值T进行收缩去噪前,主炮信号一般是连续相干的,干扰炮信号在延迟时间作用下会以串扰噪声的形式存在,所以此时分离问题可以看作噪声去除问题:

(6)

其中YXN分别表示变换域中包含噪声的混叠数据、不含噪声的原始主炮数据和串扰噪声数据.在每一次迭代中,混叠数据在Contourlet域被分解为K层,每层L个方向,用k=1, 2…, K表示分解的第k个尺度,频率越高k越大;用l=1, 2…, L表示第k个尺度的方向数.在Contourlet域,假设原始主炮数据和串扰噪声数据的各方向子带系数都服从高斯分布,根据贝叶斯最大后验概率估计理论,可以求得Wiener滤波公式(Moulin and Liu, 1999)

(7)

为第k尺度第l方向的子带中位于i, j处的系数Yi, j的去噪估计值,其中,σN2为噪声方差,σX2为原始主炮数据在第k尺度第l方向的子带系数方差.σN2σX2分别反映噪声和有效信号的大小,其比值可以表示信噪比的强弱.Wiener滤波算子的函数形式如图 1所示.Wiener滤波过程等同于将变换域中含噪系数乘以一个小于1的值,使系数实现收缩从而达到滤波目的.

图 1 Wiener阈值收缩函数形式示意图 Fig. 1 Schematic diagram of the Wiener threshold contraction function

由Wiener阈值收缩原理(式(7)),可得Wiener阈值

(8)

式(8)中σN可以采用鲁棒中位数估计法(Donoho and Johnstone, 19941995)求取

(9)

其中Fmedian为变换域最细尺度子带系数集绝对值{|Fi, j|}的中值.σX2的估计公式(李东明等,2015)为

(10)

其中NumYi, j所处子带的全部系数个数.

由于同一个子带中全部有效信号系数方差σX2的取值相同,我们记式(8)得到的阈值算子为子带一致Wiener阈值算子,在Contourlet域使用该阈值算子迭代反演分离同时震源数据的方法称为Contourlet域子带一致Wiener阈值法,下文简称为子带一致阈值法.

1.3 自适应Wiener阈值

子带一致阈值法能够取得一定的分离效果,但是由于没有考虑到串扰噪声在Contourlet域中不同尺度的分布规律以及有效信号的空间分布规律,该方法在串扰噪声的压制效果、有效信号尤其是弱信号的保留和迭代收敛速度方面存在不足.基于以上考虑,我们分别从噪声方差计算和有效信号方差计算两个方面进行改进,获得自适应Wiener阈值.

在迭代初期,由于干扰炮产生的串扰噪声很大,阈值T的选取主要考虑对串扰噪声的压制效果.此时,串扰噪声集中分布在Contourlet域的中频和高频子带中,尺度越高,噪声越多.为了加大对噪声的压制力度,不宜使用式(9)得到的全局一致噪声方差σN2,应该在高频子带中设置更大的噪声方差估计量,从而获得更小的Wiener阈值T来加大系数的收缩,提升对噪声的压制效果.因此,我们对噪声方差的计算作出如下改进:

(11)

其中,Fmax为Contourlet变换域最细尺度子带系数集绝对值{|Fi, j|}的最大值,α∈(0, 1)为衰减系数.在式(11)中,λ1为调整因子,λ2为加速因子,σN, λ1, λ22反映了不同迭代次数下不同尺度子带中噪声的待滤除水平.分解的尺度越高,频率越高,λ1的值也越大,使得σN, λ1, λ22的值也越大,从而具有尺度自适应能力;具有指数衰减特点的λ2能够提高σN, λ1, λ22逼近0的速度,进而加速Wiener阈值T趋于1,达到加快迭代收敛的目的(Gao et al., 2010Yang et al., 2013).

在迭代后期,阈值T的选取需要考虑对有效地震信号特别是弱信号的保护.由于Contourlet变换具有能量聚集性,在迭代后期噪声能量远小于连续相干信号能量的情况下,变换域的高能量区域系数主要对应有效信号.在这些系数上使用式(10)得到的子带一致σX2来估计有效信号的方差,会存在比真实值偏小的问题,导致Wiener阈值T的值偏小而加大了系数的收缩量,进而影响待分离波场中有效信号特别是弱信号的保留.因此,我们改进有效信号方差的计算方法,定义如下:

(12)

其中,Bi, j(w)表示以Yi, j为中心大小为w×w的窗口,W为窗口内系数个数,本文w取值为3.在式(12)中,有效信号方差σi, j2的计算利用了系数的局部空间位置关系,改子带一致的有效信号方差估计为加窗的局部估计,能够更加准确反映真实方差的大小,具有空间自适应的特点.

根据式(11)与式(12),改进后的尺度与空间自适应Wiener阈值算子可以表示为

(13)

在Contourlet域使用式(13)的阈值算子迭代反演分离同时震源数据的方法称为Contourlet域自适应Wiener阈值法,下文简称为自适应阈值法.

2 算例 2.1 Marmousi模型数据

Marmousi模型是石油工业界用于方法测试的一个复杂地质模型.图 2a图 2b分别为主炮和干扰炮激发采集的未混叠共检波点道集数据,共767道;图 2b数据在0~1 s随机延迟时间作用下以伪随机噪声的形式加到主炮波场中,形成如图 2c所示的混叠波场,混叠度m=2.模拟计算在Intel Xeon E7-4807、内存为128G的IBM服务器上展开,没有并行计算.Contourlet变换选择“9-7”塔式分解滤波器和“pkva”方向滤波器进行三层分解,方向数分别为4、8和8.在本次测试的自适应Wiener阈值算子中,α的取值为0.6.图 3a—b分别为图 2a数据与图 2c数据在Contourlet域中不同尺度不同方向的系数分布情况,可以看出在迭代初期,有效信号的系数集中分布在中低频部分,而串扰噪声的系数主要分布在中高频部分.图 4a图 2a数据在Contourlet域中不同尺度不同方向的系数统计分布直方图,可以看出系数多集中分布在零点附近,形态上表现为尖峰长尾,可以用高斯分布来逼近,当分解的尺度越低,得到数值较大的系数越多,进一步证明了连续相干的有效信号集中分布于低频和中频子带中;图 4b为串扰噪声数据在Contourlet域中不同尺度不同方向的系数统计分布直方图,可以看出串扰噪声分布满足高斯模型,且数值较大的系数多集中存在于中高频子带中.

图 2 Marmousi模型中共检波点道集数据 (a)原始主炮数据; (b)原始干扰炮数据; (c)混叠数据. Fig. 2 Data of Marmousi model in common receiver domain (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.
图 3 Contourlet域不同尺度不同方向分解的子带系数分布图 (a)原始主炮数据; (b)混叠数据. Fig. 3 Distribution of subband coefficients in different scales and different directions in Contourlet domain (a) Original data of the main source; (b) Blended data.
图 4 Contourlet域的不同尺度不同方向子带系数统计分布直方图 (a)原始主炮数据; (b)串扰噪声数据. Fig. 4 Statistical distribution histogram of subband coefficients in different scales and different directions in Contourlet domain (a) Original data of the main source; (b) Crosstalk noise data.

在相同的迭代反演框架下,分别用子带一致阈值法、Curvelet域指数衰减阈值法(祖绍环等,2016)和本文提出的自适应阈值法分离图 2c中的混叠数据,证明自适应阈值法的优越性.图 5a—c为主炮的分离结果,图 6a—c为干扰炮的分离结果.从分离结果可以看出,这三种方法都取得了较好的分离效果,但是自适应阈值法能够保护有效信号特别是深部弱信号,残留噪声更少.为了更好地比较这三种方法的分离效果,分别将图 2a图 5a—c数据相减、将图 2b图 6a—c数据相减得到残差,然后将残差放大10倍后得到如图 5d—f图 6d—f的放大后残差图.从残差图可以看出,子带一致阈值法没能很好地压制串扰噪声,同时也造成了有效信号的泄漏,残差大,分离结果较差;Curvelet域指数衰减阈值法对深部弱信号保留较好,但是串扰噪声强的部位不能够被很好地压制;本文提出的自适应阈值法能够有效的压制串扰噪声,更好地恢复和保护深层弱信号,分离结果好,残差小.为了更加细致的对比分离效果,我们分别提取原始主炮数据(图 2a)、混叠数据(图 2c)和分离结果(图 5a—c)中第110和300道的数据,绘制如图 7所示的单道波形对比图.从单道波形对比图中可以发现,这三种方法都能去除混叠数据中的噪声,分离的结果都与原始未混叠数据接近,但从图中放大部位可以看出本文自适应阈值法分离得到的数据与原始主炮数据更相近.图 8为从残差图 5d—f中抽取的第110和300道数据,反映了分离误差的大小,能够进一步展示图 7中的三种方法分离得到的结果与原始主炮数据的接近程度.图 8中的数值越小,证明分离结果与原始主炮数据越接近,分离效果越好.可以看出,自适应阈值法分离得到的误差比子带一致阈值法和Curvelet域指数衰减阈值方法分离得到的误差小很多,在4s以后的弱信号部位,自适应阈值法得到的分离误差几乎为0,由此证明了本文自适应阈值法具有良好地保护弱信号的能力.

图 5 不同方法分离主炮的结果和放大10倍的残差 (a)子带一致阈值法分离结果; (b) Curvelet域指数衰减阈值法分离结果; (c)本文提出的自适应阈值法分离结果; (d)、(e)和(f)分别是分离结果(a)、(b)和(c)所对应的残差(已放大10倍). Fig. 5 Deblended results of the main source from different methods and corresponding residuals with 10 times scaled up (a) Deblended result from the subband consistent threshold method; (b) from the exponential attenuation threshold method in Curvelet domain; (c) from the proposed adaptive threshold method; (d), (e) and (f) are the residuals with 10 times scaled up corresponding to the results (a), (b) and (c), respectively.
图 6 不同方法分离干扰炮的结果和放大10倍的残差 (a)子带一致阈值法分离结果; (b) Curvelet域指数衰减阈值法分离结果; (c)本文提出的自适应阈值法分离结果; (d)、(e)和(f)分别是分离结果(a)、(b)和(c)所对应的残差(已放大10倍). Fig. 6 Deblended results of the interference source from different methods and corresponding residuals with 10 times scaled up (a) Deblended result from the subband consistent threshold method; (b) from the exponential attenuation threshold method in Curvelet domain; (c) from the proposed adaptive threshold method; (d), (e) and (f) are the residuals with 10 times scaled up corresponding to the results (a), (b) and (c), respectively.
图 7 原始主炮数据、混叠数据和用三种分离方法得到的分离后主炮数据的单道波形对比图 右上角子图为红色虚线矩形框的放大图,‘A’表示子带一致阈值法,‘B’表示Curvelet域指数衰减阈值法,‘C’表示本文提出的自适应阈值法. Fig. 7 The single trace waveform comparison of the original data, blended data, and deblended data from the three deblending methods for the main source The upper right subgraph is the zoomed-in section of the red rectangle box. 'A' is for the subband consistent threshold method, 'B' for the exponential attenuation threshold method in Curvelet domain, and 'C' for the proposed adaptive threshold method. (a) The 110th trace in Fig. 2a, Fig. 2c and Figs. 5a—c; (b) The 300th trace in Fig. 2a, Fig. 2c and Figs. 5a—c.
图 8 三种分离方法得到单道分离误差对比图 蓝色虚线表示子带一致阈值法分离得到的误差,绿色实线表示用Curvelet域指数衰减阈值法分离得到的误差,红色实线表示本文自适应阈值法分离得到的误差. (a) 图 5d—f中第110道数据; (b) 图 5d—f中第300道数据. Fig. 8 A comparison of the errors for deblended single trace obtained from the three deblending methods The blue dotted line is the error for the subband consistent threshold method, the green solid line for the exponential attenuation threshold method in Curvelet domain, and the red solid line for the proposed adaptive threshold method. (a) The 110th trace in Figs. 5d—f; (b) The 300th trace in Figs. 5d—f.

为了定量分析分离效果,我们定义如下信噪比公式:

(14)

其中s表示未混叠的原始数据,表示分离结果.图 9a—c为使用这三种阈值方法分离混叠数据得到的主炮数据的分离信噪比曲线,收敛后的信噪比分别为22.57 dB、27.02 dB和30.78 dB.子带一致阈值法分离得到的信噪比较低,在2200次实现收敛,迭代次数较多,收敛慢,用时为5601.7 s,计算开销最大;Curvelet域指数衰减阈值法获得了较高的分离信噪比,迭代35次实现收敛,收敛较快,用时792.1 s,计算开销较大;本文提出的自适应阈值法有很高的分离信噪比和收敛速度,在第25次实现了收敛,用时64.5 s,用时最少.通过以上对比,进一步证明本文提出的自适应阈值法对串扰噪声有良好的压制效果,能很好地保护有效信号,用时短,收敛快,优于其他两种方法.

图 9 不同分离方法处理Marmousi模型数据的分离信噪比曲线 (a)子带一致阈值法; (b) Curvelet域指数衰减阈值法; (c)本文提出的自适应阈值法. Fig. 9 The SNR of the Marmousi model data from the different deblending methods (a) The subband consistent threshold method; (b) The exponential attenuation threshold method in Curvelet domain; (c) The proposed adaptive threshold method.
2.2 SEAM模型数据

SEAM模型包含倾斜地层、盐丘等复杂地质构造,是工业界常用的另一个测试模型.图 10a图 10b分别为主炮和干扰炮激发采集的未混叠的共检波点道集数据,共1024道,可以看到盐丘产生的反射层信号能量很强.在0~1 s随机延迟时间作用下,干扰炮波场会以伪随机噪声的形式叠加到主炮波场中,混叠度m=2,形成如图 10c的混叠波场.本次测试的分离方法是本文提出的自适应阈值法,加速因子λ2α为0.7,其他测试条件与前文相同.图 11a为主炮分离结果,收敛后的信噪比为27.88 dB,图 11b为干扰炮分离的结果,收敛后的信噪比为31.31 dB,可以看出,使用自适应阈值法分离SEAM模型的同时震源数据同样具有较高信噪比,有效信号几乎被完全恢复,位于盐丘内部的弱反射信号得到了很好的保留,同样证明了自适应阈值法能够很好地保护弱信号,具有良好的分离效果;图 11c图 10a图 11a之间的残差,图 11d图 10b图 11b之间的残差,可以看出,除了反射波能量过大的部位残留少量噪声外,大部分串扰噪声得到了很好地压制,残差较小.图 12为分离得到的信噪比曲线,通过40次迭代实现收敛,用时93.2 s,效率高.

图 10 SEAM模型数据 (a)原始主炮数据; (b)原始干扰炮数据; (c)混叠数据. Fig. 10 Data of the SEAM model (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.
图 11 本文提出的自适应阈值法分离结果和残差 (a)主炮分离结果; (b)干扰炮分离结果; (c)图 10a图 11a之间残差; (d)图 10b图 11b之间残差. Fig. 11 Deblended results and the corresponding residuals from the proposed adaptive threshold method (a) Deblended result of the main source; (b) Deblended result of the interference source; (c) Residual between Fig. 10a and Fig. 11a; (d) Residual between Fig. 10b and Fig. 11b.
图 12 本文提出的自适应阈值法分离SEAM模型数据的信噪比曲线 其中实线代表主炮的分离信噪比,虚线代表干扰炮的分离信噪比. Fig. 12 The SNR with respective to the iteration number for SEAM model data using the proposed adaptive threshold method The solid line is for the separated SNR of the main source, and the dashed line for the separated SNR of the interference source
2.3 第一个实际数据

使用Petroleum Geo-Services(PGS)公司公开的海上实际数据测试本文提出的自适应阈值法的分离效果.为了获得更加复杂的同时震源数据,我们截取原数据进行了人为混叠.实际数据的炮数和检波点数都是256,采样间隔为4 ms,混叠度m=2,加速因子λ2α为0.7,其他条件与前文相同.图 13a为原始主炮数据;图 13b为原始干扰炮数据;原始干扰炮数据在2 s内的延迟时间编码下与原始主炮数据混叠得到混叠数据,如图 13c所示.可以看出,在原始主炮数据与原始干扰炮数据中,浅部信号的能量都很强,远大于深部弱信号.

图 13 第一个实际数据 (a)原始主炮数据; (b)原始干扰炮数据; (c)混叠数据. Fig. 13 The first field data (a) Original data of the main source; (b) Original data of the interference source; (c) Blended data.

使用自适应阈值法分离图 13c中的混叠数据.图 14a为主炮分离结果,收敛后的信噪比为23.37 dB,图 14b为干扰炮分离的结果,收敛后的信噪比为21.27 dB,图 14c图 13a图 14a之间的残差,图 14d图 13b图 14b之间的残差.图 15为分离的信噪比曲线.自适应阈值法取得了较高的分离信噪比,能很好地压制串扰噪声,保护深层弱信号,由残差数据(图 14c—d)可见深部几乎没有有效信号泄漏,本次测试在第37次迭代时实现了收敛,用时89.2 s,证明了本文自适应阈值法对实际的同时震源数据同样具有优越的分离效果和高效的分离速度.

图 14 本文提出的自适应阈值法分离第一个实际数据的分离结果和残差 (a)主炮分离结果; (b)干扰炮分离结果; (c)图 13a图 14a之间残差; (d)图 13b图 14b之间残差. Fig. 14 Deblended results of the first field data from the proposed adaptive threshold method, and the corresponding residuals (a) Deblended result of the main source; (b) Deblended result of the interference source; (c) Residuals between Fig. 13a and Fig. 14a; (d) Residuals between Fig. 13b and Fig. 14b.
图 15 本文提出的自适应阈值法分离第一个实际数据的信噪比曲线 其中实线代表主炮的分离信噪比,虚线代表干扰炮的分离信噪比. Fig. 15 The SNR with respective to the iteration number for the first field data using the proposed adaptive threshold method The solid line is for the SNR of the main source, and the dashed line for the SNR of the interference source.
2.4 第二个实际数据

使用某工区实际采集的同时震源地震数据进一步验证自适应阈值法的有效性.检波点间距和炮点间距均为25 m,每道共1020个采样点,采样间隔为4 ms,混叠度m=2.在本次测试中,加速因子λ2α为0.7,其他测试条件与前文相同.图 16a为混叠的共偏移距道集,串扰噪声散布在整个道集中,深部弱有效信号被湮没.使用本文提出的自适应阈值法得到的主炮和干扰炮的分离结果分别如图 16b图 16c所示.从分离结果可见浅部强反射同相轴连续性更好,深部弱反射同相轴清晰可见,强噪声被完全压制,剖面的信噪比显著提高.本次测试进一步证明了自适应阈值法的有效性.

图 16 第二个实际数据的分离测试 (a)混叠数据; (b)主炮分离结果; (c)干扰炮分离结果. Fig. 16 The deblending test of the second field data (a) Blended data; (b) Deblended result of the main source; (c) Deblended result of the interference source.
3 讨论

在自适应阈值法中,有效信号方差σi, j2采用了局部加窗的计算方法,具有空间自适应的特点.我们用算例中Marmousi模型数据来讨论不同窗口取值大小对混叠数据分离效果的影响.窗口长度w取不同值时,使用自适应阈值法分离图 2c中的混叠数据,迭代收敛后,主炮与干扰炮的分离信噪比如图 17所示.当w取3~7时,主炮与干扰炮的分离信噪比都较大,且w取3时具有最大的分离信噪比.当w超过7时,随着窗口取值的增加,分离信噪比会持续下降.因此,w取3~7是一个合理的选择.

图 17 不同窗口长度w对本文自适应阈值法分离效果的影响 其中实线代表主炮的分离信噪比,虚线代表干扰炮的分离信噪比. Fig. 17 The effect of different window sizes w on the separation effect from the proposed adaptive threshold method The solid line is for the SNR of the main source, and the dashed line for the SNR of the interference source.

在自适应阈值法中,衰减系数α能够促进和加快迭代的收敛,我们同样使用Marmousi模型数据讨论α对分离效果的影响.当α的取值分别为0.2,0.3,0.4,0.5,0.6,0.7,0.8和0.9时,使用自适应阈值法分离图 2c中的混叠数据得到如图 18所示的主炮分离信噪比曲线.在图 18中,随着α的减小,收敛需要的迭代次数减少,收敛速度明显加快,但是当α小于0.6时,收敛后的信噪比也出现了明显的下降;当α取值在0.7~0.9之间时,收敛后的信噪比一致且较高;当α为0.9时,收敛后信噪比较高,但是收敛需要的迭代次数明显多于α为0.8时收敛需要的迭代次数.通过讨论可以发现,α越小,收敛速度越快,但是过小的α会影响分离效果,所以我们需要在分离效果和收敛速度之间做出一个平衡,而α取值为0.6~0.8是一个较为合理的选择.

图 18 衰减系数α取不同值时,本文自适应阈值法得到的主炮分离信噪比曲线 Fig. 18 The SNR curves during the separation of the main source from the proposed adaptive threshold method by taking different values of attenuation parameter α
4 结论

在同时震源数据的迭代反演分离过程中,待处理地震数据经过Contourlet变换后,有效信号和串扰噪声在不同尺度和不同方向的子带中具有不同的分布规律,因此使用子带一致Wiener阈值方法往往得不到好的分离结果.在稀疏域迭代反演框架下,本文提出了一种在Contourlet域分离同时震源数据的自适应Wiener阈值方法.在迭代初期,自适应Wiener阈值方法中的噪声方差计算具有尺度自适应能力,能够适当增大高频子带中的噪声方差估计量,相比于使用全局一致的噪声方差计算方法,更利于达到快速压制串扰噪声的目的.在迭代后期,地震数据中的有效信号在稀疏域中多表现为数值较大的能量聚集点,自适应Wiener阈值方法使用局部加窗法计算有效信号方差,比子带一致的有效信号方差能更加准确地反映实际方差的大小,能够避免因为有效信号方差计算偏小而导致弱信号泄露.实验结果证明,本文Contourlet域自适应Wiener阈值方法能够快速有效地压制串扰噪声和保护有效信号,取得了比Contourlet域子带一致Wiener阈值方法和Curvelet域指数衰减阈值方法更好的分离效果.

在估计噪声方差时,本文自适应Wiener阈值方法只考虑了噪声数据在Contourlet域中不同尺度下的分布规律,如何将同一尺度下的不同方向特征也考虑进去需要进一步研究.

References
Abma R L, Manning T, Tanis M, et al. 2010. High quality separation of simultaneous sources by sparse inversion.//80th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Akerberg P, Hampson G, Rickett J, et al. 2008. Simultaneous source separation by sparse radon transform.//78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2801-2805.
Burt P, Adelson E. 1983. The Laplacian pyramid as a compact image code. IEEE Transactions on Communications, 31(4): 532-540.
Chang S G, Yu B, Vetterli M. 2000a. Adaptive wavelet thresholding for image denoising and compression. IEEE Transactions on Image Processing, 9(9): 1532-1546.
Chang S G, Yu B, Vetterli M. 2000b. Spatially adaptive wavelet thresholding with context modeling for image denoising. IEEE Transactions on Image Processing, 9(9): 1522-1531.
Chen Y K, Fomel S, Hu J W. 2014. Iterative deblending of simultaneous-source seismic data using seislet-domain shaping regularization. Geophysics, 79(5): V179-V189.
Chen Y K, Jin Z Y, Gan S W, et al. 2015. Iterative deblending using shaping regularization with a combined PNMO-MF-FK coherency filter. Journal of Applied Geophysics, 122: 18-27.
Chen Y K. 2017. Fast dictionary learning for noise attenuation of multidimensional seismic data. Geophysical Journal International, 209(1): 21-31.
Cheng J K, Sacchi M D. 2015. Separation and reconstruction of simultaneous source data via iterative rank reduction. Geophysics, 80(4): V57-V66.
Do M N, Vetterli M. 2001. Pyramidal directional filter banks and curvelets.//Proceedings 2001 International Conference on Image Processing. Thessaloniki, Greece, Greece: IEEE, 3: 158-161.
Do M N, Vetterli M. 2002. Contourlets: a new directional multiresolution image representation.//Conference Record of the Thirty-Sixth Asilomar Conference on Signals, Systems and Computers, 2002. Pacific Grove, CA, USA, USA: IEEE.
Do M N, Vetterli M. 2003. Framing pyramids. IEEE Transactions on Signal Processing, 51(9): 2329-2342.
Donoho D L, Johnstone J M. 1994. Ideal spatial adaptation by wavelet shrinkage. Biometrika, 81(3): 425-455.
Donoho D L. 1995. De-noising by soft-thresholding. IEEE Transactions on Information Theory, 41(3): 613-627.
Donoho D L, Johnstone I M. 1995. Adapting to unknown smoothness via wavelet shrinkage. Journal of the American Statistical Association, 90(432): 1200-1224.
Doulgeris P, Bube K, Hampson G, et al. 2012. Convergence analysis of a coherency-constrained inversion for the separation of blended data. Geophysical Prospecting, 60(4): 769-781.
Eslami R, Radha H. 2003. The contourlet transform for image denoising using cycle spinning.//The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers. Pacific Grove, CA, USA, USA: IEEE, 2: 1982-1986.
Gao J J, Chen X H, Li J Y, et al. 2010. Irregular seismic data reconstruction based on exponential threshold model of POCS method. Applied Geophysics, 7(3): 229-238.
Han L G, Tan C Q, Lv Q T, et al. 2013. Separation of multi-source blended seismic acquisition data by iterative denoising. Chinese Journal of Geophysics (in Chinese), 56(7): 2402-2412. DOI:10.6038/cjg20130726
Huo S D, Luo Y, Kelamis P. 2009. Simultaneous sources separation via multidirectional vector-median filter. Geophysics, 77(4).
Kim Y, Gruzinov I, Guo M H, et al. 2009. Source separation of simultaneous source OBC data.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 51-55.
Li C B, Mosher C C, Morley L C, et al. 2013. Joint source deblending and reconstruction for seismic data.//83rd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 82-87.
Li D M, Gai M Y, Li C R, et al. 2015. Research on adaptive optics image denoising algorithm based on the Wavelet-based Contourlet transform. Laser & Optoelectronics Progress (in Chinese), 52(11): 91-96.
Li Q, Ban X G, Gong R B, et al. 2018. 2D deblending using the multi-scale shaping scheme. Journal of Applied Geophysics, 148: 44-54.
Lin T T Y, Herrmann F J. 2009. Designing simultaneous acquisitions with compressive sensing.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Liu Q, Han L G, Li H J. 2014. Synchronous interpolation and denoising in simultaneous-source data separation. Chinese Journal of Geophysics (in Chinese), 57(5): 1647-1654. DOI:10.6038/cjg20140527
Liu Z J, Wang B, Specht J, et al. 2014. Enhanced adaptive subtraction method for simultaneous source separation.//84th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 115-119.
Magnussen F. 2015. Iterative source separation using a hybrid cadzow and median filter.//85th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Moore I, Dragoset B, Ommundsen T, et al. 2008. Simultaneous source separation using dithered sources.//78th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 2806-2810.
Moore I. 2010. Simultaneous sources-Processing and applications.//80th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Moulin P, Liu J. 1999. Analysis of multiresolution image denoising schemes using generalized Gaussian and complexity priors. IEEE Transactions on Information Theory, 45(3): 909-919.
Song J W, Li P M, Qian Z P, et al. 2019. Simultaneous vibroseis data separation through sparse inversion. The Leading Edge, 38(8): 625-629.
Song J W, Li P M, Wang W C, et al. 2019. High-productivity blended acquired data separation by sparse inversion. Oil Geophysical Prospecting (in Chinese), 54(2): 268-273.
Tan C Q, Han L G, Lv Q T, et al. 2013. Separation of blended data by sparse inversion based on the reciprocity theorem. Journal of Seismic Exploration, 22(3): 209-222.
Wei Y J, Zhang P, Xu Z. 2019. Separation of 3D blending seismic data based on sparse constrained inversion. Chinese Journal of Geophysics (in Chinese), 62(10): 4000-4009. DOI:10.6038/cjg2019M0501
Wilbur J E, McDonald R J, Stack J. 2009. Contourlet detection and feature extraction for automatic target recognition.//2009 IEEE International Conference on Systems, Man and Cybernetics. San Antonio, TX, USA: IEEE, 2734-2738.
Yang P L, Gao J H, Chen W C. 2013. On analysis-based two-step interpolation methods for randomly sampled seismic data. Computers & Geosciences, 51: 449-461.
Yang Q J, Mao W J, Tang H H, et al. 2017. Deblending with weak signal preserved by dip vector-median filter.//87th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 5044-5048.
Zhang C, Olofsson B. 2012. Separating simultaneous source data using weighted tau-p transform.//82nd Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Zhang P, Mao W J. 2018. Weighted structure-enhancing constrained least-squares reverse time migration of simultaneous-source data. Chinese Journal of Geophysics (in Chinese), 61(10): 4088-4099. DOI:10.6038/cjg2018L0251
Zhang X P, Desai M D. 1998. Adaptive denoising based on SURE risk. IEEE Signal Processing Letters, 5(10): 265-267.
Zhou Y T, Li S H, Xie J Y, et al. 2017. Sparse dictionary learning for seismic noise attenuation using a fast orthogonal matching pursuit algorithm. Journal of Seismic Exploration, 26(5): 433-454.
Zu S H, Zhou H, Chen Y K, et al. 2016. Irregularly sampled multi-source data deblending based on shaping regularization. Oil Geophysical Prospecting (in Chinese), 51(2): 247-253.
Zu S H, Zhou H, Mao W J, et al. 2017. Iterative deblending of simultaneous-source data using a coherency-pass shaping operator. Geophysical Journal International, 211(1): 541-557.
Zu S H, Zhou H, Mao W J, et al. 2018. 3D deblending of simultaneous source data based on 3D multi-scale shaping operator. Journal of Applied Geophysics, 151: 274-289.
韩立国, 谭尘青, 吕庆田, 等. 2013. 基于迭代去噪的多源地震混合采集数据分离. 地球物理学报, 56(7): 2402-2412. DOI:10.6038/cjg20130726
李东明, 盖梦野, 李超然, 等. 2015. 基于小波域的Contourlet变换法的自适应光学图像去噪算法研究. 激光与光电子学进展, 52(11): 91-96.
刘强, 韩立国, 李洪建. 2014. 混采数据分离中插值与去噪的同步处理. 地球物理学报, 57(5): 1647-1654. DOI:10.6038/cjg20140527
宋家文, 李培明, 王文闯, 等. 2019. 基于稀疏反演的高效混采数据分离方法. 石油地球物理勘探, 54(2): 268-273.
魏亚杰, 张盼, 许卓. 2019. 基于稀疏约束反演的三维混采数据分离. 地球物理学报, 62(10): 4000-4009. DOI:10.6038/cjg2019M0501
张攀, 毛伟建. 2018. 同时震源数据的加权结构增强最小二乘逆时偏移. 地球物理学报, 61(10): 4088-4099. DOI:10.6038/cjg2018L0251
祖绍环, 周辉, 陈阳康, 等. 2016. 不规则采样的多震源数据整形正则化分离方法. 石油地球物理勘探, 51(2): 247-253.