地球物理学报  2017, Vol. 60 Issue (9): 3518-3538   PDF    
基于稀疏约束和多源激发的地震数据高效采集方法
王汉闯1,2, 陶春辉1,2, 陈生昌3 , 杜泳4, 丘磊1,2     
1. 国家海洋局第二海洋研究所, 杭州 310012;
2. 国家海洋局海底科学重点实验室, 杭州 310012;
3. 浙江大学地球科学学院, 杭州 310027;
4. 浙江农林大学暨阳学院, 浙江诸暨 311800
摘要:地震勘探目标日趋复杂化和精细化,"两宽一高"等采集技术获得了广泛应用,从而导致当前地震数据采集周期越来越长、成本越来越高,如何解决日益增长的勘探成本问题成为当前地震采集领域的研究热点之一.针对上述问题,本文首先开展了基于稀疏性的地震数据高效采集方法理论研究,对地震数据稀疏性基本理论、稀疏约束下随机采样及其数据重建方法进行了深入探讨,提出使用改进的分段随机采样方法灵活地进行实际地震采集测网设计;详细阐述了多源地震激发方法,对多源地震数据分离方法开展了深入研究,提出了基于小窗口中值滤波与稀疏约束联合随机去噪的多源数据分离方法,并在数据分离处理中取得了较好的效果;将上述两种地震数据采集方案有机结合,提出了1)规则多源、随机检波点(DmsRg)、2)随机多源、规则检波点(RmsDg)和3)随机多源、随机检波点(RmsRg)等三种高效采集方案及相应的数据重建方案,满足了后续常规化数据处理的要求,并讨论了多源激发对数据成像的影响.基于Marmousi模型数据的数值试验表明,本文构建的基于稀疏约束和多源激发的高效采集方法理论对于提高地震数据采集效率、降低勘探成本具有重要的应用价值,建立的数据重建方法流程可以取得和常规数据接近的成像结果.本文方法虽然在数值试验中取得了较为理想的效果,但还需要得到野外实际数据采集的进一步检验.
关键词: 勘探成本      稀疏约束      多源激发      高效采集      数据重建     
Highly efficient methods of seismic data acquisition based on sparse constraints and blended-source excitation
WANG Han-Chuang1,2, TAO Chun-Hui1,2, CHEN Sheng-Chang3, DU Yong4, QIU Lei1,2     
1. Second Institute of Oceanography, State Oceanic Administration, Hangzhou 310012, China;
2. Key Laboratory of Submarine Geosciences, State Oceanic Administration, Hangzhou 310012, China;
3. School of Earth Sciences, Zhejiang University, Hangzhou 310027, China;
4. Jiyang College, Zhejiang A & F University, Zhejiang Zhuji 311800, China
Abstract: Facing the increasingly complex and sophisticated situation of exploration targets, many new seismic acquisition technologies such as wide-azimuth, wide-bandwidth and high-density seismic data acquisition methods have been widely used, leading to the acquisition period is becoming longer and longer and the acquisition cost is becoming higher and higher than before. How to solve such a problem is a current focused issue in this field. This paper presents a study on methods and theories of highly efficient seismic data acquisition, focusing on seismic data sparse representation theory, highly efficient acquisition and reconstruction methods. We present an improved piecewise-random sampling method for seismic data sampling, and then develop a high-precise reconstruction approach for sparse acquisition. Based on the introduction of blended seismic acquisition, the blended seismic data separation is studied thoroughly. A new data separation method based on small window median filtering and sparse constraint joint random denoising is proposed to improve the practicability of the blended seismic acquisition method. We combine the above two useful efficient seismic acquisition methods and propose three highly efficient acquisition network designing methods, including 1) dense blended sources and random geophones (DmsRg) scheme, 2) random blended sources and dense geophones (RmsDg) scheme, and 3) random blended sources and random geophones (RmsRg) scheme. The corresponding reconstruction procedures are proposed which can reconstruct the seismic data on the highly efficient acquisition network to the regular and dense acquisition network. It satisfies the requirements of subsequent processing in accordance with the existing conventional means. Furthermore, the influence of simultaneous shooting on seismic data's imaging results is analyzed. We use the international standard Marmousi model to test the proposed data acquisition methods. The numerical examples show that the theoretical framework of the joint highly efficient seismic data acquisition established in this paper can improve the efficiency of the current seismic data acquisition and reduce the exploration costs. High-precise imaging result of the reconstructed data is compatible to that of conventional methods. Although the proposed methods have achieved great progress in numerical tests, more work of field seismic data acquisition is needed to validate their effectiveness in the future.
Key words: Exploration costs      Sparse constraints      Blended sources excitation      Highly efficient acquisition      Data reconstruction     
1 引言

当前地震勘探(特别是石油地震勘探)中地震数据采集的成本越来越高,这主要是因为:1) 常规地震数据采集方法是基于Fourier变换和Nyquist采样定理(Marks, 1991; Oppenheim and Schafer, 1989)的,对地震信号的采样间距(空间和时间)具有一定的要求,必须进行规则密集的信号采样;2) 随着勘探目标复杂程度的提高和后续处理解释要求的提高,地震数据需要进行“两宽一高”(宽频带、宽方位和高密度)采集.因此,地震数据高效采集方法的研究成为当前地震勘探领域的一个重要且迫切的问题.

压缩感知是一种在矩阵分析、拓扑几何、概率论、运筹学和泛函分析等学科的基础上建立起来的信息获取与处理的全新的理论体系,其基本原理是在假定信号具有稀疏性的前提下,通过特定的采样方法获取信号“精挑细选”的、极少量的采样数据,最后通过求解一个最优化问题实现原始信号的高精度重建.压缩感知最大的意义在于它突破了基于Nyquist采样定理的常规信号采样方法的局限性,实现了信号的高效采样和高分辨率重建,从而为解决信号采样中信息冗余、效率低等问题开辟了新的途径.因压缩感知利用的是信号的稀疏性,因此在本文中又称其为基于稀疏性的信号采样理论或稀疏采样理论.该理论一经提出便在医疗成像(Lustig et al., 2007)、无线通信(Mamaghanian et al., 2011)与图像超分辨重建(Park et al., 2003; Yang et al., 2010)等领域得到了高度关注,并取得了相关研究成果.在地震勘探中的应用主要集中在地震数值模拟、数据压缩、数据重建和波场非线性反演等方面(Gan et al., 2016; Hennenfent et al., 2008; Herrmann, 2010; Kim et al., 2015; Liang et al., 2014; Lin et al., 2008; Liu et al., 2015; Sternfels et al., 2015; Wang et al., 2010; Zhang et al., 2014; 唐刚, 2010; 王薇等, 2013).近年来,压缩感知也逐步被用于发展地震数据高效采集方法,主要探索如下:Milton等(2011)在常规的规则密集布置炮点的地震数据采集方式的基础上提出了随机稀疏布置炮点的数据采集方式(他们使用多维插值方法来重构数据);马坚伟(2011)的“稀疏促进地震勘探”观点中包含了利用压缩感知降低野外采集数据量的构想;陈生昌等(2015)报告了应用压缩感知理论进行地震数据高效采集方法研究的初步框架;Mosher等(2012)基于压缩感知理论提出了一种非均匀优化采样的地震数据观测系统设计方法,表明使用该方法可以获得更高的空间分辨率或者可以完成更大面积的采集任务;王汉闯(2015)王汉闯等(2016)对基于稀疏性的地震数据高效采集方法及其实施策略进行了探索性研究.以上研究成果为基于地震数据的稀疏性开展高效采集研究提供了基本思路,也为进一步开展新的高效采集方法的研究奠定了基础.

提高地震数据采集效率的另外一种途径是使用多个震源同步激发.早在20世纪70年代人们便开始研究多源激发方法(Ward et al., 1990),即使用多个(组、套)震源进行激发,当前陆地采集中常用的同步相位编码(Simultaneous Phase Encoding, SPE)(Bagaini, 2010)、远距离分隔同步扫描(Distance Separated Simultaneous Sweeping, DS3)(Bouska, 2010)以及独立同步扫描(Independent Simultaneous Sweeping, ISS) (Howe et al., 2008)等采集方式以及海上数据采集中使用的“首尾”式和“并排”式的多个震源船同步激发的采集方式(Aaron et al., 2009; Hampson et al., 2008)都属于多源方法的范畴,大大提高了地震数据采集效率.Berkhout等(2008)提出使用地震混合(blending)来表示多源激发方法,把常规的各个单震源组合起来形成一个混合源(blended sources),得到的来自多个震源的地震记录称为混合地震记录(blended data),震源混合一般采用随机延迟激发的方式.他们同时提出了多源地震数据的两种处理方法:1) 先进行混合数据的分离,得到各个单震源对应的地震数据,然后就可以按照常规的数据处理流程进行处理;2) 直接对这种混合地震数据进行偏移成像或反演处理,可以节约数据分离处理的时间.在第一种数据处理方案中,多源地震数据的分离是一个关键问题,人们发展了反问题求解方法(Ikelle, 2007)、基于稀疏性的反问题求解方法(Wason et al., 2011)、变换域滤波法(Huo et al., 2009; Mahdad et al., 2011)以及预测滤波法(Kim et al., 2009; Spitz et al., 2008)等数据分离方法来满足后续处理要求.在第二种数据处理方案中,基于震源随机时间延迟的多源地震数据的次震源干扰可以通过叠加和偏移进行压制,但在复杂地震数据中仍不能满足需求,最小二乘偏移被认为是一种有效的多源数据直接处理方法(Dai et al., 2012; Verschuur et al., 2009),在增强成像分辨率的同时可以压制串扰噪声和采集脚印(Tang et al., 2009),但过于庞大的计算量又限制了该方法的应用.

综合分析以上两种高效采集方法可知,基于稀疏性的地震数据高效采集方法是一种正在发展且有着广阔前景的采集方法,它利用了地震数据稀疏性的本质属性;多源激发方法也是近年来逐步发展起来的一种采集方式,它利用了地震波场的线性叠加原理.二者有着类似的思路,即通过一定的采集策略达到减少布置观测系统和采集数据量的目的,然后通过数据重建(稀疏恢复或者多源数据分离)得到对应常规测网的规则、密集地震数据.若把两种高效采集方法结合起来,就可以发展一种新的采集方法,可以进一步提高采集效率.本文拟在基于稀疏性的地震数据高效采集方法的基础上,在震源端把常规激发改为多源激发,从而形成一种新的高效采集方式,从1) 地震数据的稀疏性理论依据以及适合地震数据时间-空间特征的最稀疏变换,2) 基于稀疏性和多源激发的高效采集测网(测线、测点分布)设计方案和3) 高效采集数据的重建和成像策略等3个方面开展研究,最终提出了基于稀疏性和多源激发的高效采集方法,通过数值试验证明了本方法的有效性.

2 方法原理 2.1 基于稀疏性的地震数据随机采集与重建

压缩感知在地震数据高效采集方法理论的发展过程中发挥了重要作用.基于压缩感知发展地震数据高效采集方法理论的过程主要包括以下三点:1) 地震数据是否具有稀疏性?2) 根据稀疏性的先验知识来确定采样方案;3) 由采集数据必须能够准确地重构出原始数据.

2.1.1 地震数据的稀疏性

地震数据的稀疏性是利用压缩感知开展高效采集方法研究的基础.地震波场的传播满足波动方程,因而根据其解的一般表示式对地震数据的稀疏性问题进行讨论是一条合理的途径.密度均一介质、脉冲震源激发条件下的声波波动方程为

(1)

其中,r0=x0i+y0j+z0k为震源位置,u为地震波场,-f(t)δ(rr0)为震源函数.该方程的Green函数G(r, t|r0, t0)可以表示为(杜世通, 2008):

(2)

式(2) 表明,空间r0处的一个单位脉冲震源在t0时刻开始激发,所产生的球面波在t>t0时间向四周空间扩散.

待研究信号的基本特征与变换的基本波形越接近,则信号就可以取得越稀疏的变换系数(Chen et al., 1998; Starck et al., 2010).Green函数式包含了一个空间上以为基础的退化函数(即具有空间上的光滑性)和一个时间上的子波波形函数(即具有时间上的暂态性).该性质可以用于指导人们对地震波场的稀疏表示,通过选择变换基和地震波前形态相类似的一些数学变换,就可以捕捉到地震波曲线状的奇异特征,从而为地震数据提供最稀疏的表示.曲波变换(Candès et al., 2006; Candès and Donoho, 2004, 2005a, 2005b)的变换基由不同尺寸和方向的曲线状元素组成,具有与地震波前的形态非常相似曲线状的基元,近年来在地震数据的稀疏表示中占据了重要的地位(Candès and Donoho, 2004; Candès et al., 2006).

2.1.2 稀疏约束下的随机采样和重建

地震数据的采集过程可以表示为:

(3)

其中,fM为原始地震数据,dn为采集到的地震数据(nM),Rn×M是采样矩阵.设Cf的稀疏变换,即f=CHm,则式(3) 可以改写为:

(4)

其中,上标H代表Hermitian转置,mNf的稀疏变换系数,N为变换域系数的维度(若C为正交变换,则N=M;若C为冗余变换,N>M),An×N.

研究表明,对于矩阵A的任何一个子集S(记为AS,为n×|S|阶子矩阵,SA中所有可能的|S|≤k列元素的组合,|S|是S的基数,kNm中非零元素的个数),只要找到一个受限等距参数0<δk < 1使其满足:

(5)

则仍然可以对欠定问题(4) 进行求解.其中,δkA的第i列.

这表明,解m的获取依赖于AS的列非相干性.数学上证明,服从独立一致性分布的高斯随机矩阵满足此要求,即矩阵A为随机矩阵,且当稀疏度k满足

(6)

时,式(5) 总是成立的(Herrmann, 2010).

也就是说,只需要利用按照采样比n/kD·log2N获取的采样数据就可以对稀疏域中的k个非零系数进行恢复.在变换矩阵C确定的情况下,采样矩阵R也必须是随机矩阵.

虽然高斯随机采样方法适合高效采样数据的重建(Hennenfent and Herrmann, 2008),但容易造成采样点过于聚集或者分散的情况,引起相应区域采样信息的冗余或者缺失,给数据恢复处理带来一定的困难.为解决上述问题,人们发展了多种随机且均匀的采样方法,如泊松碟随机采样(Dunbar and Humphreys, 2006; Tang et al., 2009; 唐刚, 2010)、Jittered采样(Hennenfent and Herrmann, 2008)和分段采样(Wang et al., 2011)等,非常适合后续重建处理.

王汉闯等(2016)通过对比研究褶积矩阵L=AHA发展了一种改进的分段采样方法,不仅具有“蓝色噪声”频谱特征(Ostromoukhov et al., 2004),还解决了当总样点不是采样点数整数倍时出现剩余样点的问题(图 1),主要过程如下:

图 1 (a)常规分段随机采样方法,(b)改进的分段随机采样方法 Fig. 1 (a) Conventional piecewise random sampling method; (b) Improved piecewise random sampling method

1) 对所有样点分段.假设总的样点为n0,需要采样点数为ns,先对总的样点分为ns段,每个段内有int(n0/ns)个样点,剩余的样点数为nr(图 1a虚线椭圆区域),很明显nrns.里int()表示取整.

2) 段内样点数调整.从ns个段中随机选取nr个段,并把分段后剩余的nr个样点放入选好的nr个段中,调整后会有nr个段的样点数目增加1,如图 1a箭头所示.

3) 采样.从ns个段中各随机地抽取一个样点作为采样点,共ns个采样点.

该采样方法可以根据实际情况按任意采样比例(采样的道数和总道数之比)进行采样,提高了实际地震数据采集系统设计的灵活性,本文将采用此改进的分段随机采样方法开展后续研究.

通过改进的分段采样方法进行数据高效采样,采样数据必须能够实现原始数据的重建.这时需要求解下面的稀疏优化问题获得m

(7)

式中,~表示信号的估计,范数的定义为‖m0:=#{m, mi≠0},是信号稀疏性的度量,mi表示向量m的第i项.

然后,利用Lagrange乘子法把问题(7) 转化为下面的无约束优化问题,并使用范数代替范数以降低问题的求解难度,可得

(8)

其中,λ>0为拉格朗日乘子.

为了能避免矩阵求逆问题,且能够适合大规模计算,本文基于迭代阈值方法(Daubechies et al., 2004),引入加速项βn(Blumensath, 2012; Daubechies et al., 2008),提高了迭代过程的稳定性,其迭代格式如下:

(9)

其中,i为迭代次数,θ为阈值,θ0为初始阈值,并随着迭代次数按指数衰减,βn为迭代步长,保证了迭代过程的稳定性,Γimi的支撑集,Tθ为硬阈值算子,其定义如下:

(10)

通过式(9) 的迭代就可以得到满足一定精度要求的全局最优稀疏解,进而得到重建地震数据的估计值.

由此可见,在地震数据稀疏性的先验条件下,利用一定的随机采样方法进行数据采集,通过这些少量的采集数据可以重建原始数据,达到高效采集的目的.本节理论可以作为发展多种地震数据高效采集方法的基础.

2.2 多源激发采集方法理论

近年来发展起来的多源地震方法(Bagaini, 2006, 2010; Beasley, 2008; Berkhout et al., 2008; Hampson et al., 2008)可以很大程度上提高采集效率、降低勘探成本,受到了国内外的广泛关注,在陆地和海洋地震勘探中获得了大量应用.本节从多源地震方法的理论出发,研究多源地震波场传播及多源观测系统设计等问题.

2.2.1 多源激发基本理论

声波方程(1) 在频率域可以写为:

(11)

其中,rs为震源的位置,r=(x, z)为地下位置,v(r)为声波速度,u(r, ω)为位移场或压力场,S(rs, ω)为震源项.

多源是用多个震源在不同的位置上随机延迟激发的线性组合.激发时,一组震源激发完毕,然后分别移动到新的位置上进行下一组激发.设ns为常规观测系统总的震源个数,Ns为一组混合源中单源的个数,则总的多源激发数nb=ns/Ns.多源可以表示为:

(12)

其中,Sbl(rs, ω)为nb个混合激发震源,S(rs, ω)为所有ns个震源;Γnb×ns为多源混合算子,其元素Γi, k表示第i次同步激发中混合算子对第k个震源的作用,其表达式如下:

(13)

其中,fi, k为0或1,用来表示该震源是否参与同步激发(如第i次激发时,若震源k参与激发则fi, k为1,反之为0),tki为第i次激发中第k震源Sk的随机延迟时间,不参与激发则设为0.

由式(11) 和(12) 可知,多源激发条件下的波动方程可以表示为:

(14)

rg为检波器位置,由式(11) 和(14) 可知,多源Sbl激发条件下接收到地震记录Pbl(rg, ω)=ubl(rg, ω)和常规单源激发条件下的地震记录P=u(rg, ω)之间的关系可以表示为:

(15)

从式(15) 可以看出,描述常规单震源和多震源激发的地震波场的波动方程在形式上是一致的,多源地震数据是常规单源地震数据的线性叠加,叠加算子就是多源的混合算子.若把震源项改为混合震源,常规单源波动方程就可以转化为多源激发条件下的波动方程,也就可以描述相应地震波场的传播规律.

这里通过一个例子进行说明.如图 2a所示层状速度模型,整个测区(0~4000 m)设置400个震源、400个检波器,检波器排列的范围和震源排列的范围相同,炮间距和道间距都为10 m,每次使用2个震源一起激发,共可获得200个多震源地震记录(ns=400, Ns=2, nb=200).以第101个多震源(1000 m处的s1震源点和3000 m处的s2震源点)为例,两个震源的延迟激发时间分别为368 ms和789 ms,记录长度为4 s,得到的模拟地震记录如图 2b所示,是一种混合地震记录.

图 2 多源激发示例 (a)层状速度模型,从上到下层速度依次为3000、3500和4000 m·s-1;(b)含2个震源激发的混合地震记录,两个震源的延迟激发时间分别为368和789 ms. Fig. 2 Examples of blended excitation (a) Layered velocity model. Velocities from top to bottom are 3000, 3500 and 4000 m·s-1; (b) Blended seismic data with 2 simultaneous sources, whose time delays are 368 and 789 ms, respectively.
2.2.2 多源数据分离方法

多源激发得到的地震记录是一种混合记录(图 2b),必须进行分离处理才能得到各个单源记录.把式(15) 简写为如下的形式:

(16)

其中,Pblnb×1代表Pbl(rg, ω),Pns×1代表P(rg, ω),Γnb×ns为多源混合算子,满足nbns.

由多源地震波场Pbl得到各个单源地震波场P(即多源地震数据分离)的实质就是求解一个线性反问题,在数学上解非唯一,必须寻找合适的解决方案.该问题求解与Γ是否可逆以及nsnb相对大小有密切的关系.这里先把各个含噪声的地震记录的延迟时间去除,该过程称为伪解混(Mahdad et al., 2011),可以用如下方程来表示:

(17)

其中,ΓHΓ的共轭转置,它的作用是先做数据的膨胀,得到对应各个单源的地震数据,然后去除了各个地震数据的随机延迟激发时间.由此可知:伪解混结果P′不仅包含了欲得到的P,也包含了其他参与混合的震源的波场Pcross(称为交叉噪声波场),即

(18)

由于随机延迟激发时震源编码不是正交的,通过伪解混处理,交叉噪声波场Pcross在其他的道集(如共检波点道集、共中心点道集和共偏移距道集等)中并不满足线性叠加原理,交叉噪声变为随机噪声.

图 2b所示多源数据为例,经过伪解混处理的两个记录如图 3所示,处于同一反射层的同相轴顶点图 3a中的点A和图 3b中的点B处于同一时间线上(注意,消除了延迟时间后,各个记录变为3 s).这里再把P′分选到共检波点域,在原道集中没有显现出来交叉噪声波场在分选道集中就表现为随机噪声,如图 3c3d所示.

图 3 所示的含2个震源激发多震源地震数据分离图示 (a)和(b)是伪解混结果,(c)和(d)是分选到共检波点域中的结果,(e)和(f)分别是(c)和(d)中值滤波的结果,(g)和(h)分别是(c)和(d)小窗口中值滤波与稀疏约束联合去噪结果,(i)和(j)是(a)和(b)所对应的最终分离结果.各图都去除了随机延迟激发时间,时间变为3 s. Fig. 3 Illustrations of blended seismic data separation with 2 simultaneous sources (a) and (b) Pseudo-deblended results; (c) and (d) Results sorted to common receiver point domain; (e) and (f) Median filtering results of (c) and (d); (g) and (h) Median filtering and sparse constraints joint denoising results of (c) and (d); (i) and (j) Final separation results corresponding to (a) and (b). Random time delays have been removed, and the record length becomes 3 s.

在人工分选的道集中对伪解混波场进行去噪处理,任何能用来消除随机噪声的滤波方法都可应用.在人工分选的道集中进行去噪处理后,再进行抽道集就可得到消除了交叉噪声波场(这里Pcross的近似)后的伪解混波场.如果交叉噪声波场Pcross能得到完全消除,即=Pcross,则就是分离后的单震源波场.

2.2.3 基于小窗口中值滤波与稀疏约束联合随机去噪的多源数据分离方法

中值滤波是一种非线性的信号处理方法,在随机去噪方面发挥了重要作用.在一定条件下,中值滤波可以克服线性滤波处理数据(如图像)细节模糊的问题,但去噪效果依赖于滤波窗口的大小及参与中值计算的点数:窗口太小,去噪效果不好;窗口太大,又会损失太多的细节.在使用中值滤波去除地震记录孤立的随机噪声过程中,为了保证不影响到真实信号,滤波窗口不能太大,由此带来去噪不彻底的现象.图 3所示例子中,使用中值滤波对图 3c3d去噪的结果如图 3e3f所示,可以看到滤波结果中仍有部分残余噪声.

依据地震数据稀疏性原理,地震数据中的真实信号和残余的部分随机噪声在稀疏域中对应不同的系数,因而可以在稀疏域中进行分离.设f为不含噪声的地震数据,y为含噪地震数据,ε为噪声,则三者具有以下关系:

(19)

由地震数据稀疏表示理论可知,曲波变换是地震数据的一种最稀疏表示的变换,地震数据在曲波域中对应着大系数,而随机噪声和曲波的基原子具有明显不同的特征,在曲波域中对应的系数比较小.因此,可以对含噪地震数据y变换后的系数作阈值处理(去除小系数),再做曲波反变换就可以得到去噪后的结果.由此可以得到:

(20)

其中,C表示曲波变换,CH表示曲波反变换,为去噪后的地震数据,Tθ()见公式(10).

基于以上两种随机噪声去除方法,这里提出一种小窗口中值滤波与稀疏约束联合去除随机噪声的方法,即先使用小窗口中值滤波去除大部分的随机噪声,然后把中值滤波结果使用稀疏约束地震数据去噪方法进一步实现去除残余噪声的目的.该过程可以表示为:

(21)

其中,Mf()表示中值滤波算子,y经过中值滤波处理的结果,为最终去噪结果.

由于第一步中值滤波已经去除了大部分随机噪声,第二步稀疏约束地震数据去噪中取比较低的阈值;两步去噪可以在去除孤立的随机噪声的同时保持地震数据中所需要的主要信息,可以取得更好的去噪效果.最后使用小窗口中值滤波与稀疏约束联合去除随机噪声的方法对所有的共检波点道集中的随机噪声进行去噪处理,再经过抽道集即可得到最终的分离结果.

图 3的例子中,使用小窗口中值滤波与稀疏约束联合去除随机噪声的方法处理的结果如图 3g3h所示,可以看出后者的处理效果更好,孤立噪声得到了更好的压制与消除.图 3a3b两个单震源地震记录经过分离处理后的结果分别如图 3i3j所示,可以看到本文方法对多源地震数据分离取得了比较好的效果.

2.3 基于稀疏约束和多源激发的高效采集模型

本文2.1节发展了一套相对完善的地震数据高效采集方法,2.2节介绍了多源地震激发方法,并发展了一种高精度数据分离方法.本节将把上述方法理论与地震观测系统设计相结合,对高效采集的具体实施策略进行探讨.

图 4 常规规则密集采集和高效采集的震数据矩阵示意图 矩阵中黑点表示数据矩阵的元素,为某一个频率成分的复值. Fig. 4 Schematic diagram of seismic data matrix corresponding to conventional and highly efficient acquisition Black dots represent elements of data matrix, one of which is a complex value of a certain frequency component.

Berkhout指出,频率域规则地震数据(2D或3D)可以通过一个数据矩阵P来表示(Berkhout, 2012),P图 4所示,其元素Pi, j表示和检波器i、震源j对应的复值.设ng为所有的检波器数目,ns为所有的激发炮数,P的每一行代表一个共检波点道集,每一列代表一个共炮点道集,每条对角线代表一个共偏移距道集,每条反对角线代表一个共中心点道集.频率域常规单震源地震数据采集模型(Berkhout et al., 2008)可以写为:

(22)

其中,ω为圆频率,rs为震源位置,rg为检波点位置,S为震源,每一列代表一个震源排列,X(rg, rs, ω)为地下波场传播算子,G(rg, ω)代表检波器矩阵,代表一个频率分量检波器特性,每一行代表一个检波器排列,P(rg, rs, ω)为接收到的地震数据.

基于地震数据采集中的炮点和检波点之间的关系,按照稀疏性约束的要求设计随机矩阵对震源矩阵S和检波器矩阵G加以改变,就可以得到不同参数的地震数据观测系统.

2.3.1 规则多源、随机检波点(DmsRg)高效采集

按照稀疏性约束的地震数据采样矩阵设计的要求,设计检波器随机采样矩阵Θgng′×ng,其中ng′为选取的检波器数目,且ng′<ng.通过Θg作用于式(22) 中的G可得:

(23)

使用多个震源一起激发,只需要用式(12) 代替(23) 中震源项即可,此时得到的地震数据Pmg为:

(24)

(24) 式就是DmsRg高效采集的数学模型,Pmg为接收到的地震数据.

2.3.2 随机多源、规则检波点(RmsDg)高效采集

按照稀疏性约束的地震数据采样矩阵设计的要求,设计炮点随机采样矩阵为Θsns′×ns,其中ns′为按照一定的随机策略选取的炮点数目(且ns′<ns).通过Θs对式(22) 中的S加以改变,可得:

(25)

对选择的震源激发方式进行改造,则有

(26)

这里,Γ′∈nb′×ns为多震源混合算子,每一行和一个混合源相关,其元素Γkl表示震源编码(相位或者振幅编码等),nb′为多源个数,Pms为高效采集地震记录.式(26) 就是随机多源、规则检波点的高效采集方法的数学模型.

2.3.3 随机多源、随机检波点(RmsRg)高效采集

该采集方式实际上是前面两种采集方式的综合,即对激发炮数和检波器数目都进行随机采样,并在震源端使用多源激发.按照稀疏性约束的地震数据采样矩阵设计的要求,在式(26) 的基础上再使用检波器随机采样矩阵Θgng′×ng对检波器进行随机布置,可得:

(27)

(27) 式就是RmsRg高效采集的数学模型,Pmsg为高效采集地震数据.

定义DmsRg、RmsDg和RmsRg高效采集方法的效率参数μmgμmsμmsg如下:

(28)

该参数可以从理论上用来描述高效采集方法的优势,参数值越大,采集效率就越高.

分析上述过程可知,以上三种高效采集方法不仅可以提高地震数据采集效率,而且可以降低存储与传输地震数据的成本,从而可以大大降低勘探成本.根据上述高效采集方法的数学模型,本文提出使用地震数据高效采集的具体步骤如下:

1) 常规规则密集测网设计:根据勘探目标的情况和勘探工作的要求,应用常规的观测系统设计方法设计规则密集的观测测网(测线、测点和源线、源点),作为下一步进行高效采集测网设计的基准测网;

2) 高效采集测网设计:保持步骤1) 中源点的规则密集网格不动,应用稀疏约束的高效采集随机矩阵设计方法设计检波点采样矩阵和(或)炮点检波点采样矩阵(王汉闯等, 2016),在常规震源激发的基础上设计多源激发的方案(包括各个多源中同步激发震源的个数、位置以及对应的随机延迟激发时间).具体的高效采集测网包括以下三种:① DmsRg方式,保持炮点总数不变,震源使用多源激发方式使混合源的个数少于常规单炮数,使用检波点采样矩阵对检波点进行随机采样,检波器分布随机且数目大为减少;② RmsDg方式,使用震源采样矩阵对震源进行随机采样,炮点总数得到了下降,然后使用多源激发方式又使混合源的个数更加少于常规单炮总数,检波器分布规则;③ RmsRg方式,在② 设计方法的基础上,再使用检波点采样矩阵对检波点进行随机采样,震源和检波点的分布随机且都大为减少,多源激发进一步提高了效率.上述三种测网都是非规则、大间距分布的高效采集测网;

3) 执行高效采集任务:按照步骤2) 设计出的高效采集测网进行地震数据采集,得到相应类型的地震数据.

综合分析以上提出的DmsRg、RmsDg和RmsRg等三种高效采集方法,可以看到高效采集和常规地震采集方法最主要的区别就是观测测网的布置方式.常规地震采集方法的观测网格是按照Nyquist采样定理的要求而设定的规则密集、且采样间隔(空间、时间)符合采样定理的要求的网格,高效采集方法的观测网格是在常规规则密集观测网格的基础上通过随机抽样得到的大间距分布的网格,对炮点、检波点以及两者同时随机布置分别得到了不同的高效采集方法.

2.4 高效采集地震数据的重建

因为观测系统的差异,通过2.3节提出的三种采集方法获得的地震数据就具有不同的结构,如何恢复为规则密集测网上的地震数据是本节的研究内容.

2.4.1 DmsRg高效采集地震数据重建

使用DmsRg高效采集方法采集到的地震数据是一种多源激发条件下包含稀疏地震道的地震数据.因此,把高效采集地震数据重建为规则密集测网上的地震数据的处理需要经过地震道恢复和多源混合数据分离两个步骤.主要过程(图 5)如下:

图 5 DmsRg高效采集地震数据重建流程图示 Fig. 5 Flow chart of seismic data reconstruction for DmsRg high-efficiency acquisition

1) 首先基于检波点采样矩阵,对DmsRg高效采集地震数据的缺失道填充零,然后使用2.1.2节所述稀疏性约束下的地震数据重建方法对每个缺失地震道的多源混合地震数据进行稀疏约束恢复处理,得到地震数据就是常规多源地震数据;

2) 使用2.2.3节提出的多源数据分离方法对步骤1) 重建的结果进行多源地震数据分离处理,进而可以得到各个单震源所对应的地震数据,完成了地震数据的重建工作,得到了规则源点网格、规则检波点网格上的地震数据.

2.4.2 RmsDg高效采集地震数据重建

使用RmsDg高效采集方法得到的地震数据是一种随机选取的部分震源经过组合形成的多源激发条件下的数据,需要经过混合数据分离和炮数据恢复才能得到规则密集测网上的地震数据.主要处理过程(图 6)如下:

图 6 RmsDg高效采集地震数据重建流程图示 Fig. 6 Flow chart of seismic data reconstruction for RmsDg high-efficiency acquisition

1) 对多源混合地震数据进行伪解混处理,得到各炮的伪解混数据,并把缺失的炮数据填充为零;

2) 步骤1) 得到的地震数据通过抽道集转换为共检波点道集,此时,每个共检波点道集表现为随机缺失炮的特征,并且多源地震数据的交叉噪声也转化为了随机噪声;

3) 对步骤2) 所得的共检波点道集数据,使用2.1.2节所述稀疏性约束的地震数据重建方法逐个道集恢复处理,并利用2.2.3节提出的多源数据分离方法对得到的道集数据进行混合数据分离处理;

4) 把步骤3) 的结果通过抽道集转换为共炮道集数据,此时的地震数据就是规则密集测网上的地震数据.

2.4.3 RmsRg高效采集地震数据重建

使用RmsRg高效采集方法采集到的地震数据是一种对随机选择的震源经过组合形成的多震源激发条件下包含稀疏道的地震数据,这种地震数据的处理只需要在2.4.2节所述的RmsDg地震数据重建方法之前加一个稀疏道恢复处理的步骤就可以了,主要步骤(图 7)如下:

图 7 RmsRg高效采集地震数据重建流程图示 Fig. 7 Flow chart of seismic data reconstruction for RmsRg high-efficiency acquisition

1) 基于检波点采样矩阵对RmsRg高效采集地震数据缺失地震道填充零,然后使用2.1.2节所述稀疏性约束的地震数据重建方法对每个多源地震数据(含随机缺道)进行道恢复处理;

2) 步骤1) 处理得到的地震数据就是一种RmsDg地震数据,然后按照2.4.2节所述RmsDg高效采集地震数据重建方法进行后续处理.

2.5 高效采集地震数据的偏移成像 2.5.1 成像处理思路

为证明本文提出的地震数据高效采集及其重建方法的有效性,本文使用波动方程叠前偏移成像方法对高效采集地震数据及其重建地震数据的成像效果进行对比分析,处理流程如图 8所示.

图 8 高效采集地震数据偏移成像策略示意图 Fig. 8 Migration imaging strategy of seismic data from highly efficient acquisition

地震数据采集的三个过程:1) 震源激发,2) 波场传播,3) 数据接收.从上述过程看,高效采集对常规规则密集震源测网和检波点测网(或者其中之一)做了随机采样,并在震源上使用了多源激发方法,只是炮集数量以及每个炮集中道的数量有一定的减少,并未从本质上影响到地震波场的传播,因此高效采集数据和常规规则密集地震数据的偏移成像的流程是一致的.

深度方向上用下行波方程正向外推多源波场,同时用上行波方程反向外推记录波场.设W(r, rs)表示上行波算子(陈生昌等, 2001),W*(r, rg)表示下行波算子,则波场外推过程可以表示为:

(29)

其中,Sbl(r, ω)为位置r上的震源波场,Pbl(r, ω)为位置r上的检波器波场,*表示复共轭.应用时间一致性成像原理可知,位置r处的像为:

(30)

式中,Re表示取实部运算,Ibl(r)表示位置r处的成像.

由以上可知,高效采集地震数据虽然不完整,但仍然可以按照常规方法进行外推并成像,只需要按照实际情况在偏移过程中确定实际采集数据的震源和检波点的位置即可.值得注意的是,在多源高效采集地震数据中,震源和检波点(或其中之一)是按照一定的随机策略布置的,多源的个数也会根据观测系统的不同而不同,在偏移成像过程中要根据实际情况设置总的单源个数、同步激发震源个数、多源个数以及检波点个数等进行相对应的处理.

2.5.2 多震源激发对成像结果的影响

下面从理论上分析多源激发对成像结果的影响.P中第k个单源地震记录Pk可以写为:

(31)

这里用Lk表示第k个震源的传播算子,m表示地下真实模型.把式(31) 代入式(16) 可得:

(32)

其中,nb为总的多源记录数,ns为设计的单源总个数;本文稀疏布置测网条件下,实际参与混合激发的震源数为ns′(ns′≤ns由矩阵Γ的行的非零元素决定),.偏移成像算子ΦT具有如下的形式:

(33)

则偏移成像结果就是:

(34)

从式(34) 可以看出,第一项是来自各个单震源的波场的成像结果,而第二项就是多个震源之间震源波场与记录波场间不匹配而产生的交叉干扰像.从式(34) 的第二项可以看出,若Γl, iΓl, j存在比较大的差异,即震源ij距离较远或随机延迟时间差别大,则可以通过叠加和偏移很好地压制交叉干扰像.

由于偏移结果中包含了来自各个震源波场的贡献,因此从理论上说和常规单震源偏移结果是可以比拟的.野外实际生产过程中设计观测系统时尽可能把同步激发的震源距离远一些,这样就可以明显降低交叉干扰像,通过偏移成像和叠加又压制了其中部分随机干扰,直接偏移成像可以取得和常规偏移成像相当的效果.

3 地震数据高效采集数值试验

本节将基于理论模型进行数值试验,以此来检验本文提出的高效采集方法理论的可行性与有效性.

3.1 试验模型

Marmousi模型是法国石油研究院推出的国际标准模型数据(速度模型如图 9所示),该模型具有典型且复杂的断层构造,对实际勘探区非常具有代表性.本文基于此数据开展高效采集模型试验.该模型大小为nx=737, nz=750,网格间距dx=12.5 m, dz=4 m.

图 9 Marmousi速度模型 Fig. 9 Marmousi velocity model

按照2.1节所述的高效地震数据采集的具体操作步骤,首先进行常规观测测网的设计:总炮数为400,炮间距为12.5 m,第一炮位置3375 m;采用固定检波器的观测方式,第一个检波器的位置为2750 m,检波点数为500,道间距为12.5 m,排列长度为6237.5 m,观测范围如图 9中黑色线框区域.在此基础上进行后续高效采集方法的数据采集与数据重建试验.

为了表示采集数据的重建效果,这里定义第i个道集的恢复信噪比SNRi和多个道集整体恢复信噪比SNR分别为:

(35)

其中,i=1, 2, …, nn为炮集总数,fi为第i个常规道集数据,为和fi相对应的重建结果,‖‖2范数.

3.2 DmsRg高效采集试验

本节基于Marmousi模型开展DmsRg高效采集试验.用改进的分段采样方法设计DmsRg高效采集网格如下:在规则密集测网(400炮、500道)的基础上,使用2个震源同时激发,基于改进的分段采样方法随机布置250道.震源的激发工作量降低为原来的50%,检波点布置的工作量也降低为原来的50%,理论效率系数μmg=4,大大提高了数据采集的效率.

使用DmsRg高效采集方法进行采集,得到了200个混合炮记录、每个混合炮记录含随机布置250道的地震数据,图 10a展示了其中第65~69炮数据,可以看到DmsRg数据不仅包含了来自多个震源的混合波场数据,而且还存在大量的随机缺道.经过重建处理得到了对应常规规则密集测网的400炮地震数据,各炮的恢复信噪比如图 10b所示.其中,第65混合炮(即图 10a中白色线框所示的炮)的重建情况如图 11所示.从恢复结果中可以看出:1) 大部分炮记录的恢复信噪比在7 dB以上(受多源数据分离方法的限制,边界炮误差较大),含缺失道的地震数据经过处理得到了很好的重建,同相轴恢复效果良好,只有在局部大幅值区域误差稍大;2) 地震道缺失以及多源激发导致重建结果中残留部分弱交叉噪声干扰,但整体重建误差较低,说明本文的重建方法是适用于DmsRg复杂地震记录重建的.

图 10 基于图 9所示Marmousi模型的DmsRg高效采集地震数据(第65~69混合炮)(a)及所有炮数据重建结果的信噪比(b) 规则网格为400炮、500道,2个震源同时激发,稀疏网格为200混合炮、250道,整体重建信噪比为8.02 dB. Fig. 10 DmsRg seismic data based on Marmousi model shown in Fig. 9 (shot 65~69) (a) and SNR of the reconstruction results (b) The regular acquisition network is constituted by 400 sources and 500 traces, and the sparse acquisition network is constituted by 200 blended sources (the number of simultaneous sources is 2) and 250 traces. The overall recovery SNR is 8.02 dB.
图 11 图 10中白色线框所示炮数据的重建情况,重建结果对应常规地震采集的第65炮和265炮 (a)和(b)分别为重建数据,(c)和(d)为常规方法采集的数据,(e)和(f)分别为两炮的重建误差.(a)和(b)的恢复信噪比分别为7.62和9.02 dB. Fig. 11 Reconstruction result of the data shown in the white rectangle of Fig. 10, which correspond to the conventional shots 65 and 265 (a) and (b) Reconstructed data; (c) and (d) Conventional data; (e) and (f) Reconstruction errors. Recovery SNR of (a) and (b) are 7.62 and 9.02 dB, respectively.

为了进一步验证DmsRg采集方法和本文重建方法的有效性,这里对高效采集数据和恢复数据进行了偏移成像处理,结果分别如图 12a12b所示,对比常规地震数据偏移成像结果(图 12c)可以看出:1) DmsRg高效采集地震数据的直接成像结果具有比较明显的噪声,主要分布在中浅层大范围剖面上(图中红色线框中比较明显);2) 重建地震数据的成像结果中噪声得到了很大程度的压制,剖面的信噪比和光滑度都得到了明显的提升,一些微弱的偏移假象被消除.3) 高效采集数据、重建数据和常规数据的偏移成像结果在深层差别不大,主要是因为地震波场在深层有着更多的叠加,降低了随机缺失和多源激发对成像的影响.由此可见,多源激发虽然增加了地震记录的复杂度,但具有不同于基于单源的高效采集方法(王汉闯等, 2016)的优势,成像结果中的噪声更加弱且分散.但重建数据的成像结果可以达到常规规则密集数据成像的效果,因而DmsRg高效采集方法是可行的.

图 12 DmsRg高效采集试验数据成像结果 (a)高效采集数据成像结果;(b)重建数据成像结果;(c)常规数据成像结果. Fig. 12 Imaging results of DmsRg seismic data (a) Imaging result of the highly efficient acquisition seismic data; (b) Imaging result of reconstructed seismic data; (c) Imaging result of conventional seismic data.
3.3 RmsDg高效采集试验

本节基于Marmousi模型开展RmsDg高效采集试验.按照RmsDg高效采集方法设计稀疏观测系统参数如下:在400个源点中随机选择其中的200个,每2个单源组合成多源进行激发,共计100个混合炮,检波点规则布置500道,理论效率系数μms=4,这种采集方法同样提高了数据采集的效率.

使用RmsDg高效采集方法进行数据采集,得到了100混合炮、每炮含500道的高效采集地震数据.其中第5~9混合炮如图 13a所示.由于RmsDg高效采集只是在部分炮点的基础上组合成多源进行激发的,因此采集数据相比常规采集数据丢失了大量记录信息.

图 13 基于图 9所示Marmousi模型的RmsDg高效采集地震数据(第5~9炮)(a)及其重建结果的信噪比(b) 规则网格为400炮、500道,稀疏网格为100混合炮(2个震源同时激发)、500道,整体重建信噪比为7.40 dB. Fig. 13 DmsRg seismic data based on Marmousi model shown in Fig. 9 (a) and SNR of the reconstruction results (b) The regular acquisition network is constituted by 400 sources and 500 traces, and the sparse acquisition network is constituted by 100 blended sources (the number of simultaneous sources is 2) and 500 traces. The overall recovery SNR is 7.40 dB.

按照本文RmsDg高效采集地震数据的重建方法对上述采集数据进行重建处理,得到了对应常规源点和检波点网格上的规则地震数据,各个炮集的信噪比如图 13b所示.因为只有部分震源参与了激发,因此重建效果可以从激发炮和完全重建炮的恢复情况进行分析:1) 激发炮以第5混合炮(图 13a中白色线框区域)为例,经过重建处理得到了对应的第5和211的常规单炮记录(图 14A_a和(A_b),信噪比为8.56和9.19 dB;2)图 14B为没有任何采样信息、完全重建得到的第223炮记录,信噪比为8.45 dB,噪声压制的也不错.从重建结果来看,内部区域重建较好,而突变点(大幅值区域)误差稍大;3) 重建数据和常规直接采集数据具有很高的一致性,证明了RmsDg高效采集方法和本文恢复方法的可行性.

图 14 RmsDg高效采集地震数据重建结果 (A)图 13中白色线框区域所示炮的重建结果,对应常规地震采集的第5炮和211炮,其中(A_a)和(A_b)分别为重建数据,(A_c)和(A_d)为常规方法采集的数据,(A_e)和(A_f)分别为两炮的重建误差.(A_a)和(A_b)的恢复信噪比分别为8.56和9.19 dB.(B)缺失炮的重建结果,对应常规采集的第223炮,其中(B_a)为重建数据,(B_b)为常规方法采集的数据,(B_c)为重建误差.该炮的恢复信噪比为8.45 dB. Fig. 14 RmsDg seismic data and its reconstruction results (A) Reconstructed results of the data shown in the white rectangle of Fig. 13 which correspond to shots 5 and 211. (A_a) and (A_b) are the reconstructed results (Recovery SNR: 8.56 and 9.19 dB), (A_c) and (A_d) are the conventional seismic data, (A_e) and (A_f) are the reconstructed errors. (B) Reconstructed results of the blank shot, which corresponds to shot 223. (B_a) is the reconstructed result, (B_b) is the conventional seismic data, (B_c) is the reconstruction error. The recovery SNR of (B_a) is 8.45 dB.

本文再对高效采集数据和恢复数据进行偏移成像处理,结果分别如图 15a15b所示,通过和图 15c所示的常规地震数据偏移成像结果比较可以看出:1) RmsDg高效采集地震数据的直接成像剖面会受到一些噪声的污染,这种噪声也分布在中浅层整体剖面上(图中红色线框中比较明显),和DmsRg数据成像结果类似;而经过处理得到的地震数据的成像剖面的信噪比得到了提升,接近于常规地震数据成像结果.2) 采样数据和恢复数据在深层都可以得到比较好的成像效果.由此可以证明,RmsDg高效采集方法可以提高采集效率、降低勘探成本.

图 15 RmsDg高效采集试验数据成像结果 (a)高效采集数据成像结果;(b)重建数据成像结果;(c)常规数据成像结果. Fig. 15 Imaging results of RmsDg seismic data (a) Imaging result of high efficient acquisition seismic data; (b) Imaging result of reconstructed seismic data; (c) Imaging result of conventional seismic data.
3.4 RmsRg高效采集试验

本节基于Marmousi模型开展RmsRg高效采集试验.按照RmsRg高效采集方法设计稀疏观测系统参数如下:基于改进的分段采样方法,在400个源点中随机选择其中的200个,每2个单源组合成多源进行激发,共计100个混合炮,在500个常规规则密集检波点中随机布置250道.震源激发和检波器布置的工作量各为原来的50%,而多源激发方法又使激发工作量降低50%,理论效率系数μmsg=8.从原理上看,RmsRg高效采集方法和RmsDg高效采集方法震源布置和激发方法一致,但布置了更少的检波器,其工作效率是上述3种高效采集方法中最高的,采集数据量也是最少的.

在上述观测系统的基础上,按照RmsRg高效采集的步骤进行采集,得到了100混合炮、每炮含250道的高效采集地震数据.其中第5~9混合炮如图 16a所示,可以看出RmsRg高效采集地震数据是一种缺失地震道的多源记录数据.

图 16 基于图 9所示Marmousi模型的RmsRg高效采集地震数据(第5~9混合炮)(a)及各炮数据的重建信噪比曲线(b) 规则网格为400炮、500道,稀疏网格为100炮(2个震源同时激发)、250道,炮采样率是50%,整体重建信噪比为6.12 dB. Fig. 16 RmsRg seismic data based on Marmousi model shown in Fig. 9 (a) and SNR of the reconstruction results (b) The regular acquisition network is constituted by 400 sources and 500 traces, and the sparse acquisition network is constituted by 100 blended sources (the number of simultaneous sources is 2) and 250 traces. The overall recovery SNR is 6.12 dB.

对高效采集数据进行重建处理,各炮的重建信噪比曲线如图 16b所示,由于本方法是对布置的部分炮点组合成多源进行激发的,因此经过重建得到的各炮记录的信噪比也会有差异,表现出一定的震荡.对比图 13b所示的RmsDg高效采集数据的恢复信噪比可以看出,由于RmsRg高效采集地震数据的数据量更少,从而导致各炮的恢复信噪比降低1.0~2.0 dB.以第5混合炮(图 16中白色线框区域)为例,经过重建处理得到了对应的第5和211的常规单炮记录(图 17A_aA_b),信噪比为6.86和7.13 dB,可以看到缺失地震道信息获得了不错的重建,整体剖面的连续性比较好,重建误差主要集中在大振幅同相轴周围.图 17B为没有任何采样信息、完全重建得到的第223炮记录,信噪比为6.58 dB,从重建结果来看,仍然是内部区域重建较好,而大振幅同相轴附近误差较大.从本次采集和重建结果中可以看出,RmsRg高效采集具有效率上的优势以及实际可操作性.

图 17 RmsRg高效采集地震数据重建结果 (A)图 16中白色线框区域所示炮的重建结果,对应常规地震采集的第5炮和211炮,其中(A_a)和(A_b)分别为重建数据,(A_c)和(A_d)为常规方法采集的数据,(A_e)和(A_f)分别为两炮的重建误差.(A_a)和(A_b)的恢复信噪比分别为6.86和7.13 dB.(B)缺失炮的重建结果,对应常规采集的第223炮,其中(B_a)为重建数据,(B_b)为常规方法采集的数据,(B_c)为重建误差.该炮的恢复信噪比为6.58 dB. Fig. 17 RmsRg seismic data and its reconstruction results (A) Reconstructed results of the data shown in the white rectangle of Fig. 16 which correspond to shots 5 and 211. (A_a) and (A_b) are the reconstructed results (Recovery SNR: 6.86 and 7.13 dB), (A_c) and (A_d) are the conventional seismic data, (A_e) and (A_f) are the reconstructed errors. (B) Reconstructed results of the blank shot, which corresponds to shot 223. (B_a) is the reconstructed result, (B_b) is the conventional seismic data, (B_c) is the reconstructed errors. The recovery SNR of (B_a) is 6.58 dB.

和前面两个试验一样,这里也对高效采集数据和恢复数据进行了偏移成像处理.结果见图 18.可以看到:1) RmsRg高效采集地震数据的直接成像结果的整体剖面相比常规偏移成像结果(图 18c)具有分布在中浅层剖面上、比DmsRg和RmsDg采集数据成像结果的噪声更大的干扰(图中红色线框中比较明显),这主要是由采集数据量的减少引起的;2) 重建地震数据的成像结果中噪声得到了很大程度的压制,提高了成像的信噪比,更接近于常规成像结果,表明了RmsRg高效采集方法的可行性.从成像结果对比中也可以看出,高效采集地震数据量不足是造成成像噪声的主要原因,而本文恢复重建方法则可以很好地补充缺失的炮和道数据,即使重建数据仍与常规单炮数据存在一定的误差,也能够在后续偏移中通过叠加改善成像质量,得到更接近常规地震成像结果的剖面.

图 18 RmsRg高效采集试验数据成像结果 (a)高效采集数据成像结果;(b)重建数据成像结果;(c)常规数据成像结果. Fig. 18 Imaging results of RmsRg seismic data (a) Imaging result of highly efficient acquisition seismic data; (b) Imaging result of reconstructed seismic data; (c) Imaging result of conventional seismic data.
4 结论

针对当前地震勘探中存在的低采集效率、高成本问题,基于地震数据稀疏性原理的高效采集方法获得了广泛研究和快速发展,而且地震多源激发方法在加快地震数据采集效率上也发挥了重要作用.本文将地震数据的稀疏性和多源激发地震方法相结合,对地震数据高效采集新方法进行了积极的探索研究,取得了一些有意义的成果.主要结论如下:

(1) 基于稀疏性的地震数据高效采集方法突破了以Nyquist采样定理为基础的常规地震数据采集方案的限制,使用随机采样方法布置大尺度随机分布的观测系统可以大大加快数据采集效率;

(2) 多源激发方法使用多个震源同时激发的策略,改变了地震波产生和数据采集的物理路径,其波场传播过程仍然满足波动方程.根据波场的线性叠加原理可以实现混合地震数据的分离,从而获得与常规单震源相对应的地震数据.本文发展的基于小窗口中值滤波与稀疏约束联合随机去噪的地震数据分离方法具有良好的效果;

(3) 本文将基于稀疏性的地震数据高效采集方法和多源地震激发方法有机结合,提出了一套更为高效的地震数据采集方案(包括三种高效采集策略),并实现了混合地震数据的高精度重建.数值计算表明,本文高效采集方案及重建方法可以满足后续数据处理要求,具有理论可行性和方法有效性;

(4) 高效采集地震数据的直接成像剖面上产生了一些噪声和假象,主要分布在剖面的中浅层;重建数据的成像剖面质量获得了较大提升,接近于常规地震数据成像结果.另外,采集数据和重建数据在深层都可以得到比较好的成像效果.

综上所述,本文构建的基于稀疏约束和多源激发的地震数据高效采集理论体系对于提高采集效率、降低勘探成本具有重要的应用价值.本文方法虽然在数值试验中取得了较为理想的效果,但还需要野外实际数据采集的进一步检验,这也是本文未来努力的方向.

参考文献
Aaron P, Van Borselen R, Fromyr E. 2009. Simultaneous sources:a controlled experiment on different source configurations.//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 56-60.
Bagaini C. 2006. Overview of simultaneous vibroseis acquisition methods.//76th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 70-74.
Bagaini C. 2010. Acquisition and processing of simultaneous vibroseis data. Geophysical Prospecting, 58(1): 81-100. DOI:10.1111/gpr.2010.58.issue-1
Beasley C J. 2008. A new look at marine simultaneous sources. The Leading Edge, 27(7): 914-917. DOI:10.1190/1.2954033
Berkhout A J. 2012. Seismic Migration:Imaging of Acoustic Energy by Wave Field Extrapolation. Vol.12. Elsevier.
Berkhout A J G, Blacquière G, Verschuur E. 2008. From simultaneous shooting to blended acquisition.//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2831-2838.
Blumensath T. 2012. Accelerated iterative hard thresholding. Signal Processing, 92(3): 752-756. DOI:10.1016/j.sigpro.2011.09.017
Bouska J. 2010. Distance separated simultaneous sweeping, for fast, clean, vibroseis acquisition. Geophysical Prospecting, 58(1): 123-153. DOI:10.1111/gpr.2010.58.issue-1
Candès E J, Donoho D L. 2004. New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities. Communications on Pure and Applied Mathematics, 57(2): 219-266. DOI:10.1002/cpa.v57:2
Candès E J, Donoho D L. 2005a. Continuous curvelet transform:Ⅰ. Resolution of the wavefront set. Applied and Computational Harmonic Analysis, 19(2): 162-197. DOI:10.1016/j.acha.2005.02.003
Candès E J, Donoho D L. 2005b. Continuous curvelet transform:Ⅱ. Discretization and frames. Applied and Computational Harmonic Analysis, 19(2): 198-222. DOI:10.1016/j.acha.2005.02.004
Candès E, Demanet L, Donoho D, et al. 2006. Fast discrete curvelet transforms. Multiscale Modeling & Simulation, 5(3): 861-899.
Chen S C, Cao J Z, Ma Z T. 2001. Prestack depth migration method based on quasi-linear Born approximation. Chinese Journal of Geophysics (in Chinese), 44(5): 704-710. DOI:10.3321/j.issn:0001-5733.2001.05.014
Chen S C, Chen G X, Wang H C. 2015. The preliminary study on high efficient acquisition of geophysical data with sparsity constraints. Geophysical Prospecting for Petroleum (in Chinese), 54(1): 24-35.
Chen S S, Donoho D L, Saunders M A. 1998. Atomic decomposition by basis pursuit. SIAM Journal on Scientific Computing, 20(1): 33-61. DOI:10.1137/S1064827596304010
Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration. Geophysical Prospecting, 60(4): 681-695. DOI:10.1111/gpr.2012.60.issue-4
Daubechies I, Defrise M, De Mol C. 2004. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics, 57(11): 1413-1457. DOI:10.1002/(ISSN)1097-0312
Daubechies I, Fornasier M, Loris I. 2008. Accelerated projected gradient method for linear inverse problems with sparsity constraints. Journal of Fourier Analysis and Applications, 14(5-6): 764-792. DOI:10.1007/s00041-008-9039-8
Du S T. 2008. Seismic Wave Dynamics:Theory and Method (in Chinese). Dongying: China Petroleum University Press.
Dunbar D, Humphreys G. 2006. A spatial data structure for fast Poisson-disk sample generation.//Proceedings of the ACM SIGGRAPH 2006. New York, NY, USA:ACM.
Gan S W, Wang S D, Chen Y K, et al. 2016. Compressive sensing for seismic data reconstruction via fast projection onto convex sets based on seislet transform. Journal of Applied Geophysics, 130: 194-208. DOI:10.1016/j.jappgeo.2016.03.033
Hampson G, Stefani J, Herkenhoff F. 2008. Acquisition using simultaneous sources. The Leading Edge, 27(7): 918-923. DOI:10.1190/1.2954034
Hennenfent G, Herrmann F J. 2008. Simply denoise:wavefield reconstruction via jittered undersampling. Geophysics, 73(3): V19-V28. DOI:10.1190/1.2841038
Herrmann F J. 2010. Randomized sampling and sparsity:Getting more information from fewer samples. Geophysics, 75(6): WB173-WB187. DOI:10.1190/1.3506147
Howe D, Foster M, Allen T, et al. 2008. Independent simultaneous sweeping a method to increase the productivity of land seismic crews.//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2826-2830.
Huo S D, Luo Y, Kelamis P. 2009. Simultaneous sources separation via multi-directional vector-median filter.//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 31-35.
Ikelle L. 2007. Coding and decoding:seismic data modeling acquisition and processing.//77th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 66-70.
Kim B, Jeong S, Byun J. 2015. Trace interpolation for irregularly sampled seismic data using curvelet-transform-based projection onto convex sets algorithm in the frequency-wavenumber domain. Journal of Applied Geophysics, 118: 1-14. DOI:10.1016/j.jappgeo.2015.04.007
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.
Liang J W, Ma J W, Zhang X Q. 2014. Seismic data restoration via data-driven tight frame. Geophysics, 79(3): V65-V74. DOI:10.1190/geo2013-0252.1
Lin T T, Lebed E, Erlangga Y A, et al. 2008. Interpolating solutions of the Helmholtz equation with compressed sensing.//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts.
Liu W, Cao S Y, Li G F, et al. 2015. Reconstruction of seismic data with missing traces based on local random sampling and curvelet transform. Journal of Applied Geophysics, 115: 129-139. DOI:10.1016/j.jappgeo.2015.02.009
Lustig M, Donoho D, Pauly J M. 2007. Sparse MRI:The application of compressed sensing for rapid MR imaging. Magnetic Resonance in Medicine, 58(6): 1182-1195. DOI:10.1002/(ISSN)1522-2594
Ma J W. 2011. Sparsity-promoting seismic exploration.//Proceedings of the 27th Annual Meeting of China Geophysical Society (in Chinese). Changsha:China Geophysical Society.
Mahdad A, Doulgeris P, Blacquiere G. 2011. Separation of blended data by iterative estimation and subtraction of blending interference noise. Geophysics, 76(3): Q9-Q17. DOI:10.1190/1.3556597
Mamaghanian H, Khaled N, Atienza D, et al. 2011. Compressed sensing for real-time energy-efficient ECG compression on wireless body sensor nodes. IEEE Transactions on Biomedical Engineering, 58(9): 2456-2466. DOI:10.1109/TBME.2011.2156795
Marks R J. 1991. Introduction to Shannon Sampling and Interpolation Theory. New York: Springer-Verlag.
Milton A, Trickett S, Burroughs L. 2011. Reducing acquisition costs with random sampling and multidimensional interpolation.//81th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 52-56.
Mosher C, Li C B, Kaplan S, et al. 2012. Non-uniform optimal sampling for simultaneous source survey design.//84th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts.
Oppenheim A V, Schafer R W. 1989. Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice-Hall.
Ostromoukhov V, Donohue C, Jodoin P M. 2004. Fast hierarchical importance sampling with blue noise properties. ACM Transactions on Graphics, 23(3): 488-495. DOI:10.1145/1015706
Park S C, Park M K, Kang M G. 2003. Super-resolution image reconstruction:a technical overview. IEEE Signal Processing Magazine, 20(3): 21-36. DOI:10.1109/MSP.2003.1203207
Spitz S, Hampson G, Pica A. 2008. Simultaneous source separation:a prediction-subtraction approach.//78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2811-2815.
Starck J L, Murtagh F, Fadili J M. 2010. Sparse Image and Signal Processing:Wavelets, Curvelets, Morphological Diversity. Cambridge: Cambridge University Press.
Sternfels R, Viguier G, Gondoin R, et al. 2015. Multidimensional simultaneous random plus erratic noise attenuation and interpolation for seismic data by joint low-rank and sparse inversion. Geophysics, 80(6): WD129-WD141. DOI:10.1190/geo2015-0066.1
Tang G, Shahidi R, Herrmann F J, et al. 2009. Higher dimensional blue-noise sampling schemes for curvelet-based seismic data recovery.//79th Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts, 191-195.
Tang G. 2010. Seismic data reconstruction and denoising based on compressive sensing and sparse representation[Ph. D. thesis] (in Chinese). Beijing: Tsinghua University.
Tang Y X, Biondi B. 2009. Least-squares migration/inversion of blended data.//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2859-2863.
Verschuur D J, Berkhout A J. 2009. Target-oriented least-squares imaging of blended data.//79th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2889-2893.
Wang D L, Tong Z F, Tang C, et al. 2010. An iterative curvelet thresholding algorithm for seismic random noise attenuation. Applied Geophysics, 7(4): 315-324. DOI:10.1007/s11770-010-0259-8
Wang H C. 2015. Theoretical research on methods of high efficient seismic data acquisition[Ph. D. thesis] (in Chinese). Hangzhou: Zhejiang University.
Wang H C, Tao C H, Chen S C, et al. 2016. Study on highly efficient seismic data acquisition method and theory based on sparsity constraint. Chinese Journal of Geophysics (in Chinese), 59(11): 4246-4265. DOI:10.6038/cjg20161126
Wang W, Han B, Tang J P. 2013. Regularization method with sparsity constraints for seismic waveform inversion. Chinese Journal of Geophysics (in Chinese), 56(1): 289-297. DOI:10.6038/cjg20130130
Wang Y F, Cao J J, Yang C C. 2011. Recovery of seismic wavefields based on compressive sensing by an l1-norm constrained trust region method and the piecewise random subsampling. Geophysical Journal International, 187(1): 199-213. DOI:10.1111/gji.2011.187.issue-1
Ward R M, Brune R H, Ross A, et al. 1990. Phase encoding of vibroseis signals for simultaneous multisource acquisition.//60th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts.
Wason H, Herrmann F J, Lin T T. 2011. Sparsity-promoting recovery from simultaneous data:a compressive sensing approach.//81th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 6-10.
Yang J C, Wright J, Huang T S, et al. 2010. Image super-resolution via sparse representation. IEEE Transactions on Image Processing, 19(11): 2861-2873. DOI:10.1109/TIP.2010.2050625
Zhang H, Chen X H, Li H X. 2014. 3D seismic data reconstruction based on complex-valued curvelet transform in frequency domain. Journal of Applied Geophysics, 113: 64-73.
陈生昌, 曹景忠, 马在田. 2001. 基于拟线性Born近似的叠前深度偏移方法. 地球物理学报, 44(5): 704–710. DOI:10.3321/j.issn:0001-5733.2001.05.014
陈生昌, 陈国新, 王汉闯. 2015. 稀疏性约束的地球物理数据高效采集方法初步研究. 石油物探, 54(1): 24–35.
杜世通. 2008. 地震波动力学理论与方法. 东营: 中国石油大学出版社.
马坚伟. 2011. 稀疏促进地震勘探. //中国地球物理学会第二十七届年会论文集. 长沙: 中国地球物理学会.
唐刚. 2010. 基于压缩感知和稀疏表示的地震数据重建与去噪[博士论文]. 北京: 清华大学.
王汉闯. 2015. 地震数据高效采集方法理论研究[博士论文]. 杭州: 浙江大学.
王汉闯, 陶春辉, 陈生昌, 等. 2016. 基于稀疏约束的地震数据高效采集方法理论研究. 地球物理学报, 59(11): 4246–4265. DOI:10.6038/cjg20161126
王薇, 韩波, 唐锦萍. 2013. 地震波形反演的稀疏约束正则化方法. 地球物理学报, 56(1): 289–297. DOI:10.6038/cjg20130130