2. 中国地质科学院, 北京 100037;
3. 中国石化石油物探技术研究院, 南京 210014
2. Chinese Academy of Geological Sciences, Beijing 100037, China;
3. Sinopec Geophysical Research Institute, Nanjing 210014, China
地震数据采集往往是在经济效益和采集质量之间寻求平衡,传统采集在相邻震源之间设置足够的激发间隔,以避免接收记录相互混叠,其代价是地震测量的参数(炮距、道距、方位角、偏移距分布等)都不一定能够达到最佳.与地震勘探作业链的其它环节相比,地震采集方法技术特别是采集思路的变化是缓慢的.多源地震同时激发采集(simultaneous acquisition)是地震数据记录方面的重大变革,多个不同空间位置的震源以一定编码方式同时激发,获得波场相互干涉的混叠炮记录.从时间域的间断激发、逐炮接收,到同时激发、连续接收,这是地震采集思路的重要突破,极大加快了数据记录速度,进而提升采集效率和成像质量.该采集技术在国外早已提出,并在近几年迅速发展:Silverman[1]提出可控震源的同时激发采集;Beasley等[2]引入了脉冲型震源的同时激发采集;Bagaini[3]对比讨论了各种不同可控震源同时激发采集方法.多源地震混合采集(blended acquisition)是在此基础上进一步发展起来的:Vaage[4]提出随机或者伪随机延迟激发震源的海上电缆采集;Berkhout[5]引入混合采集,即对震源项进行随机线性编码,使得不同震源按照一定的随机时间延迟激发,并将混合震源激发与混合检波器接收联合,提出双重混合的激发采集系统[6];Blacquiere[7]进一步论述了混合激发采集相对于传统采集方法的双重照明优势.国内学者佘德平[8]对多震源混合炮记录的数值模拟进行了研究;陈生昌等[9-10]对多震源数据的合成及超炮道集的偏移成像技术进行探讨和研究.
多源地震混合采集技术实现的关键在于炮分离(deblending),即将混合采集记录中相互混叠的多源波场记录分离,得到传统采集的单炮记录.利用最小二乘算法得到的伪分离记录中,混叠噪声在共炮点道集以外的其它时域呈现为随机分布,因此去噪滤波成为了一种思路:Moore等[11]和Akerberg等[12]在伪分离道集的共偏移距道集进行滤波去噪;Huo等[13]采用全方位中值滤波器在其共中心点道集压制混叠噪声;Doulgeris等[14]和Mahdad等[15]进一步设计了迭代算法,在共检波点道集提取混叠噪声,优化分离结果;Doulgeris等[16]进一步对比分析了不同滤波方法的处理效果;我们从滤波[17]和反演[18]两方面对炮分离做了一些研究.目前主要限制条件仍然在于混合度,即多源混叠数据中的震源数量,由于混叠噪声随着混合度的提高成倍增加,导致炮分离质量锐减;另外目前算法往往还需较高的迭代次数,计算成本较高.
本文针对高混合度的混叠炮记录,研究了一种多时域迭代去噪的炮分离方法,该方法主要运用了多级中值滤波技术和Curvelet阈值迭代去噪算法.前者由于其简单有效和自适应的特性,在近年来得到了迅猛发展,并被广泛应用于地震资料的去噪处理中[19-22];后者由Herrmann[23],Hennenfent[24]等引入到压制地震资料随机噪声领域,并有着出色的表现.周怀来等[25]和刘磊等[26]也对此进行了研究及应用.本文将这两种方法应用在炮分离中:在伪分离记录的共偏移距道集采用多级中值滤波消除混叠噪声,再将滤波结果转换到共检波点道集以Curvelet阈值迭代去噪法进一步压制残留噪声,最终返回共炮点道集得到分离结果;在迭代过程中需要逐步降低中值滤波的时窗长度,不断提取有效信号,将其恢复到炮分离结果中.实际资料处理结果证明,本方法能够明显提升炮分离质量和计算效率.
2 混叠数据 2.1 常规采集Berkhout[27]提出地震数据可以很方便地用数据矩阵P来描述:设计一个含N个检波点的2D地震测线,在整个观测过程中检波器位置不动.炮点放置在第一个检波器处,一炮之后炮点向前移动一个炮检距长度,最后得到N个炮记录.将采集得到的炮记录作傅里叶变换到频率域,然后按照每一个频率分量重排,重排后的数据存放在一个矩阵中(图 1).
在这个数据矩阵中,各列代表单频率的共炮点道集,各行代表单频率的共接收点道集,主对角线存放的是共偏移距道集,次对角线存放的是共中心点道集.我们可以对每一个频率分量都构建这样的数据矩阵,很多的地震数据处理中,都可以对每个频率独立进行处理,因而这种矩阵表示在地震数据处理中作用很大.模型试算采用150个震源激发,炮间距10m,150个检波器接收,道间距10m,炮点与检波点重合,采样间隔4 ms,常规逐炮采集后抽取不同时域的地震记录如图 2.
Berkhout[5]对混合激发采集技术做了详细阐述,不同空间位置的多个震源按照随机线性编码方式激发构成时域混叠的炮记录,由下式表达:
(1) |
Pbl表示混叠炮集记录中某一频率的共频道集,P表示常规采集的单炮记录相应频率共频道集,Γ表示混合算子,其列向量Γn与混合炮集中相应列的混合炮记录对应,它所包含的元素表示震源编码项。由于采用线性编码,Γmn=emn-jωτ,τmn定义了该混合炮集第n个混合炮记录中第m个单炮的延迟激发时间.在混合采集中,不同震源按照一定的随机时间延迟激发.Berkhout还定义了混合度(blendingfactor):
(2) |
它是混合采集的核心要素,实际上表示混合炮记录中所包含的震源数量,野外实际采集时该比值理论上越高越好:如果给定震源总数,那么采集效率会随混合度增加而成倍提高;如果给定采集时间,那么所记录到的数据量就会随之成倍增加,进而提高成像质量.然而这将导致不同震源信号在时域中的混叠程度加剧,炮分离难度倍增.以上述模型为例,将等距间隔500m的三个震源混合激发,总共得到50个混合度为3的混合炮记录.抽取不同时域记录如图 3所示.
单从数学上看炮分离可以看作求解(1)式,但由于混合炮Pbl个数少于单炮P的个数,因此(1)式是欠定的,只能求其最小平方解:
(3) |
由于采用线性编码方式,因此其结果相当于共轭转置.进而得到伪分离记录:
(4) |
P′即为求得的伪分离记录,其不同时域道集如图 4所示:可见,伪分离将混合炮记录按照其混合度进行复制(在该模型中即复制三份),并按照各单炮对应的延迟激发时间分别将其归位,但其他震源信号仍然完整保留,即为“混叠噪声”(blendingnoise)或“串扰噪声”(crosstalk).
混合激发采集技术实现的关键在于炮分离,也就是将相互混叠的各单炮记录分离出来,即可进行常规后续处理.该采集模型混合度为3,所以混叠噪声能量约为有效信号的2倍.由图 4可见,由于对震源采用随机线性相位编码激发,该噪声仅在共炮点道集连续分布,在其它时域(共偏移距道集、共检波点道集)均为随机分布,而有效信号在任何时域都是连续的.只要能够在某时域将混叠噪声消除,即可实现炮分离.然而问题在于,待分离信号往往淹没在噪声之中,去噪难度较大;而且混叠噪声实际上是其它震源激发产生的有效信号,与一般随机噪声(斑点噪声、椒盐噪声)有着本质的区别,许多常规去噪方法难以适用.本文运用了多级中值滤波与Curvelet阈值迭代去噪算法,在不同时域根据噪声特点采用相应的去噪方法,并通过迭代运算优化炮分离结果.
3.1 基本方法中值滤波是地震资料去噪处理中的常用工具,由于其简单有效,到目前为止,该方法已经发展成为一项比较成熟的技术.在此基础上改进的多级中值滤波器是针对图像处理提出的一种最具代表性的细节保护中值滤波器之一,这种滤波器通过一组能较好地“匹配”图像细节的基本结构---基本子窗口有效区分信号结构和噪声,从而达到保护细节结构和消除噪声的目的,而其滤波作用还会随着滤波时窗长度的增加而增加,因此被广泛用于随机噪声的压制[19-22].经过试验发现,大步长多级中值滤波对于高混合度采集数据中的混叠噪声亦具有很好的压制作用.图 5a即为在共偏移距道集11点多级中值滤波的结果.可见大部分混叠噪声能够被去除.但由于混叠程度较高,仍有部分残留噪声,因此还需要做进一步的处理.
考虑到混叠噪声在共检波点道集亦为随机分布,因此将上一步去噪结果转换到该时域,这样残留噪声不仅再次呈现随机分布特征,且经过上一步处理后,其能量已经大幅衰减,表现为一般随机噪声,进一步采用Curvelet阈值迭代去噪算法对其进行压制.该算法基本思想是将压制地震随机噪声的问题转换为L1范数最优化反演问题[23-24],其基本原理是把欲处理的数据变换到Curvelet域中,在Curvelet域中设定适当的阈值进行限制,最后将处理后的数据变换到时空域.在Curvelet域中同相轴有效信息的系数和随机噪声的系数仅有较小的重叠,这样就提供一个很好的信噪分离方法.阈值的大小可由数据的平均值确定.使用该算法将残留的混合噪声去除(图 5b),然后转换回共炮点道集即得到分离结果(图 5c).
3.2 迭代算法虽然大步长的多级中值滤波强大的平滑作用将大量的混叠噪声消除殆尽,但是有效信号也被严重破坏,大部分的深部弱反射信号和多次波都丢失.这是由于该滤波方法存在严重不足即是引起相对滤波窗口而言较为“细小”的信号细节结构的破坏和丢失.研究表明,长度较小的窗口能够较好地保护信号的细节信息,但却不能有效地滤除噪声;而长度较大的窗口能更好地抑制噪声,同时却严重地损失重要信息.在初次滤波时,为了尽可能将成倍于有效信号的混叠噪声去除,势必选择较大的滑动窗口,但大量有效信号也会不可避免的丢失或者被破坏.为了将其恢复,我们进一步将分离结果按照混合算子混合后再伪分离(图 6a),计算其与伪分离采集记录P′之间的差值:
(5) |
由图 6b可见,仍有大量细节构造信息存在于该残差中.抽取该残差数据的共偏移距道集再次进行中值滤波,此时应减小滤波长度,以便提取更多有效信号.图 7a所示为7点多级中值滤波结果,将滤波结果R′i+1加入到上一次迭代的炮分离结果:
(6) |
再转换至共检波点道集进行Curvelet阈值迭代去噪,消除残留噪声(图 7b),最后得到本次计算后的炮分离结果Pi+1(图 7c).可见,损失的有效信号得到了补偿,分离结果被进一步优化.
随着迭代次数增加,炮分离结果将被不断优化.迭代的收敛取决于分离结果与目标函数P′之间的差值,将迭代运算后得到的分离结果进行伪分离,计算其与目标数据之间的信噪比:
(7) |
式中,下标rms代表均方根,分离质量越高,该信噪比就越高.在进行炮分离之前,应首先根据实际采集的混合度和数据资料的品质确定期望值,如在某次迭代后未能达到事先设定的期望值,则将两者之间的差值Ri+1再次代入到下一步运算;如果已经达到或者该信噪比不再提高则终止运算.在本文模型试算中将期望值设定为15dB,而经过计算第二次迭代之后该信噪比约为10dB(图 8),即认为分离结果仍可继续优化,因此进行了第三次迭代运算,结果如图 9所示.可见,更多的细节信号被提取,原本被混叠噪声完全淹没的深部弱反射信号和多次波都被较好地分离出来.此时,分离结果的伪分离记录与目标数据的信噪比超过了15dB.由图 10可见,在差剖面中几乎没有残留的有效信号,因此终止迭代.
1.计算伪分离炮记录P′=PblΓH,并根据混合度和资料品质设置信噪比期望值、中值滤波初始时窗长度及其迭代步长,令i=0,Pi=0为初始炮分离结果;
2.i=i+1;
3.计算剩余混叠信号Ri=P′-Pi-1ΓΓH;
4.将剩余混叠信号选排在共偏移距道集,按照初始时窗长度和步长计算得到本次滤波的时窗长度进行中值滤波,消除混叠噪声,提取有效信号R′i;
5.将滤波结果与上次炮分离结果相加:P′i=Pi-1 +R′i;
6.将P′i转换到共检波点道集进行曲波阈值迭代去噪,压制残留噪声;
7.转换回共炮点道集得到分离结果Pi;
8.计算炮分离结果的伪分离记录PiΓΓH及其与目标数据P′的信噪比,如果达到期望值则跳出循环,输出Pi为最终炮分离结果,否则回到第二步.
4 实际资料处理实际资料处理时采用了某海上区块的实际地震资料:其采样间隔4ms,道间距和炮间距均为25m,观测时间为4s,常规采集得到150个单炮记录(图 11a),对其数值模拟合成总共15个混合度为10的混叠炮记录(图 11b).按照下式计算每一次迭代运算得到的分离结果中有效信号与混叠噪声的信噪比,以此判定炮分离的质量:
(8) |
由图 11b可见,伪分离炮记录中混叠噪声的能量远远超出有效信号,尤其多次波和深部弱反射信号几乎湮没在了噪声之中,其信噪比仅为-2.79dB.采用本文迭代去噪算法进行炮分离处理如图 11所示:第一次去噪虽将大量的混叠噪声消除殆尽,但有效信号也被严重破坏,大部分的细节构造都丢失了;在第二次处理过程中降低了中值滤波的时窗长度,部分丢失的信号被恢复,与此同时采用Curvelet阈值迭代去噪算法在共检波点道集控制噪声水平,有效压制了残留噪声;第三次去噪处理后,大部分细节构造都显现出来,分离结果被进一步优化;经过四次迭代之后,炮分离结果(图 11c)就已经与实际单炮道集十分接近,其信噪比达到了18.57db,信噪比提升约21.36db,由差剖面(图 11d)可见几乎没有有效信号残留.
5 讨论将Mahdad等[15]、Doulgeris等[16]研究结果与本方法处理结果对比如表 1,分析可见:前者方法旨在预测混叠噪声,将其从伪分离记录中不断剥离,在两炮混合激发的情况下,混叠噪声与有效信号水平相当,此时能够得到较好的分离结果;但当混合度增大时,混叠噪声随之成倍增加,预测难度陡增,即便对于三炮混合的情况,分离质量的下降也十分明显,而且其庞大的计算成本也不容忽视.本文旨在提取有效信号,首先使用大步长的多级中值滤波将混叠噪声尽可能平滑滤除,然后在迭代过程中降低滤波时窗长度,从残差数据中不断恢复细节构造信息,并在共检波点道集采用Curvelet阈值迭代算法压制误入噪声,优化分离结果.从分离质量和迭代次数来看,本方法都取得了更好的效果,能够用于分离高混合度的混叠数据.
多源地震混合采集推进了地震采集技术的进步,带来了采集理念、地震资料数据质和量的深刻变化,是未来地震采集方法发展的必由之路.尽管对于混叠数据直接进行成像处理可能是未来的发展方向,而近年来也涌现出成功的例子,但就目前而言,能够得到高质量叠前数据的炮分离技术仍然是主流方向.相对于反演炮分离,基于滤波去噪的炮分离方法更加快速实用.
混合度是混合采集的核心要素,在实际采集中,该比值越高越好,其代价是混叠噪声随之倍增,加大去噪分离难度.本文重点研究高混合度混叠数据的分离,提出了一种多时域联合迭代去噪的方法.与前人方法预测混叠噪声,将其从混叠数据中剥离不同,本方法旨在预测有效信号.对于高度混叠的混合炮记录,混叠噪声的预测要比有效信号的提取更加困难,因此通过后者实现炮分离可能是更好的选择.从实际资料处理结果来看,本方法不仅能够得到高质量的分离结果,而且大幅提升了计算效率.
[1] | Silverman D. Method of three dimensional seismic prospecting. U. S. Patent, 4159463, 1979. |
[2] | Beasley C J. A new look at marine simultaneous sources. The Leading Edge , 2008, 27(7): 914-917. DOI:10.1190/1.2954033 |
[3] | Bagaini C. Overview of simultaneous vibroseis acquisition methods. 76th Annual International Meeting, SEG, Expanded Abstracts, 2006: 70-74. |
[4] | Vaage S. Method and system for acquiring marine seismic data using multiple sources. U.S. Patent 6906981B2, 2002. |
[5] | Berkhout A J, Blacquière G, Verschuur D J. Changing the mindset in seismic data acquisition. The Leading Edge , 2008, 27(7): 924-938. DOI:10.1190/1.2954035 |
[6] | Berkhout A J G, Blacquiere, Verschuur D J. The concept of double blending: Combining incoherent shooting with incoherent sensing. Geophysics , 2009, 74(4): A59-A62. DOI:10.1190/1.3141895 |
[7] | Blacquière G, Berkhout G, Verschuur E. Double illumination in blended acquisition. 81th Annual International Meeting, SEG, Expanded Abstracts, 2011: 11-15. |
[8] | 佘德平. 多震源地震正演数值模拟技术. 地球物理学进展 , 2012, 27(4): 1533–1540. She D P. A numerical modeling of simultaneous sources. Progress in Geophys. (in Chinese) , 2012, 27(4): 1533-1540. |
[9] | 陈生昌, 马在田. 广义地震数据合成及其偏移成像. 地球物理学报 , 2006, 49(4): 1144–1149. Chen S C, Ma Z T. Generalized synthesis of seismic data and its migration. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1144-1149. |
[10] | 陈生昌, 王汉闯, 陈林. 三维VSP数据高效偏移成像的超道集方法. 地球物理学报 , 2012, 55(1): 232–237. Chen S C, Wang H C, Chen L. A high efficient super-gather migration method for 3D VSP data. Chinese J. Geophys. (in Chinese) , 2012, 55(1): 232-237. |
[11] | Moore I, Dragoset W, Ommundsen T, et al. Simultaneous source separation using dithered sources. 78th Annual International Meeting, SEG, Expanded Abstracts, 2008: 2806-2809. |
[12] | Akerberg P, Hampson G, Richett H, et al. Simultaneous source separation by sparse radon transform. 78th Annual International Meeting, SEG, Expanded Abstracts, 2008: 2801-2804. |
[13] | Huo S D, Luo Y, Kelamis P, et al. Simultaneous sources separation via multi-directional vector-median filter. 79th Annual International Meeting, SEG, Expanded Abstracts, 2009: 31-35. |
[14] | Doulgeris P, Mahdad A, Blacquiere G. Separation of blended impulsive sources using an iterative approach. 72nd Annual International Meeting, EAGE, Extended Abstracts, 2010: B004. |
[15] | Mahdad A, Doulgeris P, Blacquiere G. Separation of blended data by iterative estimation and subtraction of blending interference noise. Geophysics , 2011, 76(3): Q9-Q17. DOI:10.1190/1.3556597 |
[16] | Doulgeris P, Mahdad A, Blacquière G. Iterative separation of blended marine data: discussion on the coherence-pass filter. 81th Annual International Meeting, SEG, Expanded Abstracts, 2011: 26-31. |
[17] | Tan C Q, Han L G, Zhang Y H, et al. Separation of blended data by iterative denoising. 74nd Annual International Meeting, EAGE, Extended Abstracts, 2012: A045. |
[18] | Tan C, Han Q, Lü Y, et al. Separation of blended data by sparse inversion based on reciprocity theorem. Journal of Seismic Exploration, 2013, (in press). |
[19] | 刘财, 李红星, 陶春辉, 等. 模糊嵌套多级中值滤波方法及其在地震数据处理中的应用. 地球物理学报 , 2007, 50(5): 1534–1542. Liu C, Li H X, Tao C H, et al. A new fuzzy nesting multilevel median filter and its application to seismic data processing. Chinese J. Geophys. (in Chinese) , 2007, 50(5): 1534-1542. |
[20] | 刘洋, 王典, 刘财, 等. 局部相关加权中值滤波技术及其在叠后随机噪声衰减中的应用. 地球物理学报 , 2011, 54(2): 358–367. Liu Y, Wang D, Liu C, et al. Weighted median filter based on local correlation and its application to poststack random noise attenuation. Chinese J. Geophys. (in Chinese) , 2011, 54(2): 358-367. |
[21] | 王旭松, 杨长春. 对地震图像进行保边滤波的非线性各向异性扩散算法. 地球物理学进展 , 2006, 21(2): 452–457. Wang X S, Yang C C. An edge-preserving smoothing algorithm of seismic image using nonlinear anisotropic diffusion equation. Progress in Geophys. (in Chinese) , 2006, 21(2): 452-457. |
[22] | 张尔华, 王伟, 高静怀, 等. 非线性各向异性扩散滤波器用于三维地震资料噪声衰减与结构特征增强. 地球物理学进展 , 2010, 25(3): 866–870. Zhang E H, Wang W, Gao J H, et al. Non-linear anisotropic diffusion filtering for 3D seismic noise removal and structure enhancement. Progress in Geophys. (in Chinese) , 2010, 25(3): 866-870. |
[23] | Hennenfent G, Herrmann F J. Simply denoise: Wavefield reconstruction via jittered undersampling. Geophysics , 2008, 73(3): V19-V28. DOI:10.1190/1.2841038 |
[24] | Herrmann F J, Hennenfent G. Non-parametric seismic data recovery with curvelet frames. Geophysical Journal International , 2008, 173(1): 233-248. DOI:10.1111/gji.2008.173.issue-1 |
[25] | 周怀来, 李枚, 郑文锋, 等. 基于Curvelet变换的地震资料信噪分离技术. 地球物理学进展 , 2009, 24(2): 620–625. Zhou H L, Li M, Zheng W F, et al. An approach to signal noise separation of seismic data based on Curvelet transform. Progress in Geophys. (in Chinese) , 2009, 24(2): 620-625. |
[26] | 刘磊, 刘振, 张军华. 曲波阈值法地震弱信号识别及去噪方法研究. 地球物理学进展 , 2011, 26(4): 1415–1422. Liu L, Liu Z, Zhang J H. Denoising and detecting seismic weak signal based on curvelet thresholding method. Progress in Geophys. (in Chinese) , 2011, 26(4): 1415-1422. |
[27] | Berkhout A J. Seismic Migration: Imaging of Acoustic Energy by Wave Field Extrapolation. New York: Elsevier Scientific Pub. Com., 1982 . |