2014, Vol. 57
2. 中国石油化工股份有限公司石油物探技术研究院, 南京 210014
2. Sinopec Geophysical Research Institute, Nanjing 210014, China
采集效率与处理效率的提高是地震采集技术一直追求的目标.传统观测系统的设计受一定的放炮时间间隔门槛值的制约，施工周期都较长；为了缩短采样周期，在观测系统设计的时候总是尽可能地减少激发炮数，从而导致激发震源点不足、资料品质差.近年来国外兴起的多震源激发技术(Bagaini, 2006; Beasley et al., 1998; Moerig et al., 2002)则可以克服这个缺陷，大大地提高观测效率， 缩短采样周期，从而降低勘探成本，并已经在近几年地震数据采集中逐步得到应用(Aaron et al., 2009; Beasley, 2008; Berkhout et al., 2009; Berkhout et al., 2008; Blacquiere et al., 2009; Hampson et al., 2008)，在地震勘探中表现出了高效、低成本等巨大的优越性.多震源地震方法在海洋与陆地高效地震数据采集中都开始得到应用，所使用的震源主要为可控震源(Bagaini, 2006; Moerig et al, 2002).陆上多震源采集如独立同步扫描方法(Howe et al., 2009; Howe et al., 2008)(ISS)、远距离同步扫描方法(Bouska, 2010)(DS3)等一般用于三维大面积地震数据采集，多台震源车间隔一定的距离独立工作；为了减少彼此之间的干扰，各个震源之间一般间隔较远，极大地提高了生产率并一定程度上压制了震源之间的干扰.相比陆地勘探，海洋多震源地震方法(Beasley, 2008; Hampson et al, 2008; Ting and Zhao, 2009)应用更加广泛，震源主要使用气枪震源，放置震源的船按照一定的排列进行随机激发，放置拖缆检波器束的船同时接收来自各个震源的波场数据.这种工作方式极大地提高了工作效率，并可以进行宽方位角地震数据采集.所谓多震源就是在多个不同位置上激发震源，可以是同步激发也可以是非同步激发(如线性时间延迟、随机激发等)，因此得到的地震记录是一种由多个震源产生的混合地震记录，如何对这种多震源地震数据进行处理是亟待研究的问题.
对于多震源地震数据，若可以准确地分离为常规单震源地震数据，就可以按照现有的方法进行成像处理.多个震源激发的情况下，对某一个震源的波场来讲，来自其他震源的波场可以认为是噪音.因此可以把多震源地震数据的分离问题看作一个去噪问题来研究，通过几何滤波(Beasley, 2008; Beasley et al, 1998)、变换域滤波(Akerberg et al., 2008)、多方向矢量中值滤波(Huo et al., 2009)等方法可以达到分离多震源地震数据的目的；也可以把多震源数据的采集过程看成一个数学过程，把多震源地震数据的分离问题转化为一个反演问题来研究，采取反演(Ikelle, 2007)与迭代反演的方法(Mahdad et al., 2011)求取单震源地震数据；多次扫描(含一次扫描)多震源地震数据的分离问题(Wang et al., 2013)也得到了深入讨论.多震源地震数据分离后就得到了各个震源的地震数据，然后就可以按照常规单震源地震偏移成像方法进行处理.
在目前的单震源地震勘探方法技术中，波动方程叠前深度偏移在实际数据处理中获得了越来越多的认同，共炮点道集数据叠前深度偏移成像可以获得高质量的偏移成像结果，但因边界处理而增加计算量等原因致使它的计算效率很低.该方法的计算效率还与炮道集数有关，炮道集数越多，计算效率越低(Berkhout, 2008).针对波动方程共炮点道集偏移成像的上述特点，人们提出了一系列改进方法，如平面波偏移成像方法(Berkhout, 1992; Chen et al., 2012; Rietveld et al., 1992; Zhang et al., 1999; 陈树民等, 2001)、相位编码(Romero, et al., 2000; Sun, et al., 2002)方法与小波束方法(Wu and Chen, 2001)等. 陈生昌根据地震数据的线性叠加原理提出 了地震数据广义合成方法(陈生昌和马在田, 2006)， 又把合成函数选取为随机函数，发展为广义随机合成偏移成像方法(陈生昌等, 2012).按照这种方法，人们可以根据不同的地质情况和要求通过人工合成，得到各种不同的人工合成震源和地震数据道集，如平面波震源和数据道集、局部平面波震源和道集以及面向目标的人工合成震源和数据道集，对这些人工合成的震源与地震数据可以直接进行偏移成像处理，大大提高了地震数据偏移成像的计算效率.这些都为发展更高效、更准确、更具针对性的波动方程叠前深度偏移成像方法技术奠定了基础.与广义合成震源类似，多震源就是在多个不同位置上激发震源，是一种按照一定格式编码的多个震源，所形成的波场是多个震源波场叠加的混合波场；不同的是广义合成震源与广义合成地震数据是通过人工手段对常规单炮地震数据所做的处理，目的是为了提高室内数据处理效率，而多震源是按照一定的采集系统要求设计的，多震源地震数据是在新的采集系统下直接采集得到的.由此可见，多震源数据的偏移成像可以采取与广义合成地震数据相似的处理方法.
本文一方面把多震源地震数据分离问题作为一个数学反演问题来研究，通过多震源地震数据分离方法把多震源地震数据分离为常规单震源地震数据，然后按照常规单程波偏移成像方法进行偏移成像处理；另一方面在陈生昌提出的广义随机合成道集偏移成像方法的基础上，结合Berkhout等提出的Blended Source和Blended Data及其偏移 成像的概念(Berkhout, 2008; Berkhout et al., 2009)， 把多震源看作经过编码的广义震源，多震源地震数据也就是这种广义震源激发条件下的地震数据，提出了多震源地震数据直接偏移成像方法.把本文提出的两种处理多震源数据的成像方法应用于数值模拟的多震源地震数据中，并与常规单震源地震方法成像结果进行对比，结果表明本文方法是正确的、有效的.
二维密度均一的频率域声波方程可以表示为
这里，v(x,z)为位置 r (x,z)处的声波速度， u (x,z,ω)为位置 r 处的位移场或压力场， r s(xs,zs)为震源位置， S ( r s,ω)为震源项.
多震源是用多个震源在不同的位置上随机延迟激发的组合(即经过编码的广义震源)，表达式可以写为
这里， Γ ( r s,ω)为地震波场混合矩阵算子：
其中t( r s)是位置 r s上各个震源的随机延迟激发时间，α( r s)为与炮点坐标有关的权系数.为了保持多震源中各个震源函数的振幅一致性，一般令权系数α( r s)=1，有
由式(1)和(2)可知，多震源激发条件下的波动方程可以表示为
设 r d(xd,zd)为检波器位置，由式(1)和(5)可知，多震源 S bl( r s,ω)激发条件下检波器接收到地震记录 P bl= u bl(xd,zd,ω)和常规单震源激发条件下的地震记录 P = u (xd,zd,ω)之间的关系就可以表示为
可以看出，多震源地震数据是常规单震源地震数据的线性叠加，叠加算子就是多震源的混合矩阵算子.
同时，由式(1)和(5)也可以看出，描述常规单震源和多震源激发的地震波场的波动方程在形式上是一致的，区别主要是两者右端的震源项.把震源项改为多震源，方程就可以描述多震源激发条件下的地震波场的传播规律.
地震波场数值模拟是研究复杂介质中地震波传播规律的有效途径之一，这里在前一小节多震源地震波场传播理论的基础上进行多震源波场的数值模拟.多震源数值模拟的观测系统在常规观测系统基础上重新设计，如图1所示，设观测排列上总的震源个数为nshot，把整个震源排列分为N部分，每部分的震源个数为M(则nshot=N×M)；每个部分的中第i(i=1,2,…,M)个震源一起组成一组一起激 发的多震源(即震源1、M+1、2M+1、…、(N－1)M+1 等N个震源组成一个多震源，震源2、M+2、2M+2、…、(N－1)M+2等N个震源组成一个多震源…，共有M对组合).检波器的分布范围固定，即在测区范围固定好检波器，然后移动炮点进行激发而不再移动检波点，第一个检波点与第一震源点的距离为xl.这样，在所有M个多震源激发观测之后就可以得到M个多震源地震数据，也就完成了整个测区的观测.与常规单震源地震采集系统相比，多震源地震采集系统一次激发多个炮点，明显加快了采集效率.
由式(2)—(3)及图 1所示观测系统可以看出， 影响多震源的因素主要有3个：(1)多震源中激发震源的个数；(2)多震源中各个激发震源位置；(3)多震源中各个震源的随机延迟激发时间，设置合适的多震源参数对多震源模拟波场的特征有着重要影响.多震源地震波场数值模拟主要是改变震源函数，即使用多个震源同步激发，然后按照波动方程进行传播，最后得到多震源地震数据.
如图2a所示层状速度模型，从上到下层速度依次为3000、3500、4000 m/s，在图 1所示的观测系统中，整个测区(0~4000 m)设置400个震源、400个检波器，检波器排列的范围和震源排列的范围相同(xl=0)，炮间距和道间距都为10 m，令N=2,M=200，即每次使用两个震源一起激发，共可获得200个多震源地震记录.以第101个多震源(1000 m处的s1震源点和3000 m处的s2震源点)为例，两个震源的延迟激发时间分别为368 ms和789 ms，记录时间点为4001，长度为4 s，时间采样率为1 ms，得到的模拟地震记录如图 2b所示.可以看出，多震 源地震记录是来自两个单震源激发的混合地震记 录，同相轴相互交织在一起，形成一种复杂的地震记录.
 | 图1 多震源地震观测系统.整个观测区域分成N部分，每个部分有M个震源，每个部分第i个一起组成 一个多震源(i=1,2,…,M)，共M个多震源.xl为第一个震源与第一个检波点之间的距离，检波点位置固定不变 Fig.1 Multi-source seismic acquisition system. Partition the whole acquisition region to N subintervals, each of which contains M sources. The i-th source of each subinterval (i=1,2,…,M) form a blended sources, and M blended sources have been formed in total by this way. xl is the distance from the 1st source to the 1st receiver, and all the receivers are fixed in the acquisition system |
 | 图2 层状速度模型及其所对应的含两个震源激发的地震记录 (a)层状速度模型，从上到下层速度依次为3000 m/s、3500 m/s、4000 m/s;(b)含两个震源激发的 混合地震记录，两个震源的延迟激发时间分别为368 ms和789 ms. Fig.2 The layered velocity model and one of the blended shot gathers with 2 simultaneous sources (a) The layered velocity model with the velocity of the three layers are 3000 m/s, 3500 m/s and 4000 m/s from the top to the bottom; (b) One of the blended shot gathers with 2 simultaneous sources whose random time delays are 368 ms and 789ms respectively. |
多震源地震波场中，对于欲得到的那个震源波场而言，其他震源产生的波场可以看作是噪音，若去除了这些噪音也就实现了多震源地震数据的分离.下面结合2.2节模拟得到的200个含两个震源的多震源地震记录(其中一炮如图2b所示)为例说明多震源地震数据的分离方法.
把式(6)写为离散形式，则多震源地震数据可以表示为
式中，k=1,2,…,N，每个单震源波场 P k随机延迟激发得到的多震源地震波场 P bl是通过混合矩阵 Γ k对各个单震源波场 P k做时移然后再叠加.
由于地震波场满足线性叠加原理，多震源波场中的交叉噪音在其本来的观测道集中是显现不出来的.考虑到每个震源波场都具有一定的随机延迟时间，可以反时移消除 P bl中施加在各个单震源波场上的随机延迟时间，这一步可通过混合矩阵 Γ k预分离处理，即有
其中， Γ Hk是 Γ k的共轭转置，从上式可知：预分离结果 P Hk不仅包含了要求的 P k也包含了其他震源波场 P i(i≠k)(交叉噪音波场)，但 P Hk相对于 P k消除了随机延迟激发时间.如图 2所示模型试验中，200个含两个震源激发的地震记录数据经过预分离得到了400个单震源地震记录，每个记录都消除了延迟时间，处于同一反射层的同相轴顶点(图3a中的点A和图3b中的点B)处于同一时间线上.
由于多个震源激发的波场包含在采集到的多震源地震数据中，且随机延迟激发时震源编码不是正交的，通过预分离处理的交叉噪音波场 P i(i≠k)在其他的人工分选道集(如共检波点道集、共中心点道集和共偏移距道集等)中并不满足线性叠 加原理，就明显地表现出来了(Mahdad et al, 2011)， 随机延迟激发方式产生的交叉噪音变为随机噪音.这里把 P Hk分选到共检波点域，即预分离得到的400个地震记录分选到共检波点域，对某一欲分离数据来说则来自其他震源的波场都作为随机噪音的形式出现在道集中，如图3c、图 3d所示；
在人工分选的道集中对预分离波场进行去噪处理，任何能用来消除随机噪音的滤波方法都可以应用.在人工分选的道集中进行去噪处理后，再进行抽道集就可得到消除了交叉噪音波场 P i(i≠k)后的
预分离波场 P′ k= P Hk- P i(i≠k).如果 P Hk中交叉噪音波场 P i(i≠k)能得到完全消除，则 P′ k
就是待求的第k号震源的单震源波场 P k，即得到了分离波场 P .这里用中值滤波进行随机噪音去除，对于图 2b所示的多震源记录，分离之后得到如图 3e、图 3f所示的两个单震源地震记录.
 | 图3 图2b所示的含两个震源激发多震源地震数据分离图示 (a)和(b)是预分离结果;(c)和(d)是分选到共检波点域中的结果;(e)和(f)是 最终分离结果.图中各图都去除了随机延迟激发时间，时间变为3 s. Fig.3 The separation of the blended shot gather with 2 simultaneous sources that shown in Fig.2b (a) and (b) are the pseudo-debelending results, and that in common receiver domain are shown in (c) and (d), (e) and (f) are the final separation results. The random time delays of all the seismic shot gathers are eliminated, and the record length become 3 s. |
在多震源地震波场传播理论、数值模拟及数据分离的基础上，针对多震源地震数据的特点，本文提 出了两种多震源地震数据的偏移成像方法：分离后再偏移成像及直接偏移成像处理.
完成了多震源地震数据的分离处理之后，就可以按照常规偏移成像方法进行处理.设震源深度为z0，检波器深度为zd，深度方向上用下行波方程正向外推多震源波场，即利用波场算子(陈生昌等, 2001)W(zn,z0)将震源波场从z0传播到zn；同时用上行波方程反向外推多震源记录波场，即利用W*(zn,zd)将检波器波场从zd传播到zn，也就获得了深度zn上的震源波场 S (x,zn,ω)和检波器波场 P (x,zn,ω),这个过程可以表示为
其中*表示复共轭，应用时间一致性成像条件提取偏移成像结果，则深度zn处的像可以由如下公式给出：
式中，x表示横向位置，zn表示深度，ω表示频率，Re表示取实部运算， I (x,zn) 表示(x,zn)处的像.
多震源激发条件下，多震源地震波场满足波动方程(5)，因此可采用与常规单震源激发地震偏移成像相似的技术方案，差别在于震源波场进行传播时，震源项是多个单震源的组合，记录波场反向传播时， 使用的是混合地震记录，即来自各个震源的波场记录.
设多个震源的深度为z0，并令式(2)的深度为z0，即 r s=(x,z0)，据此就可以得到多震源地面激发的波场 S bl，记为 S bl(x,z0,ω).这个过程可以写为
其中， Γ (x,z0,ω)=α(x,z0,ω)ejωt(x,z0)，表示地面处多震源波场的混合算子.把式(11)表示的混合震源取代方程(9)的第一个式子中的震源 S (x,z0,ω)进行震源波场传播，即
由此就可以得到各个深度层zn上的多震源波场 S bl(x,zn,ω).
设检波器深度为zd，则记录波场 P bl可以写为 P bl(x,zd,ω).用 P bl(x,zd,ω)取代方程(9)中第二个式子的记录波场 P (x,zd,ω)进行检波器波场的反向传播，即使用单程波外推算子的共轭W*(zn,zd)将检波器波场从zd传播到zn，可得深度zn上的检波器波场 P bl(x,zn,ω)：
得到了深度zn上的震源波场 S bl(x,zn,ω)和检波器波场 P bl(x,zn,ω)之后，根据时间一致性成像条件，深度zn处的像可以用如下公式给出：
式中，x表示横向位置，ω表示频率，*表示复共轭，Re表示取实部运算， I bl(x,zn) 表示(x,zn)处的像.经过深度外推就可以得到各个深度的成像结果，最后输出整个成像剖面.
由以上可知，多震源成像和常规单震源成像类似，只要分别在多个震源加上随机激发时间的影响，即通过混合矩阵算子把多个位置上激发的单震源组合为多震源，然后通过波场传播得到地下各个区域的多震源波场，结合多震源地震记录波场，应用时间一致性成像原理可以得到最终的成像结果.
下面从理论上分析多震源激发对成像结果的影响.在式(7)中， P k可以写为
这里Lk是第k个震源对应的传播模拟算子.把上式代入式(7)可得
其中， m 表示地下真实模型， Φ =∑ N k=1 Γ k L k.偏移成像算子 Φ T就可以定义如下的形式：
则偏移成像结果就是
从式(18)就可以看出，第一项是来自各个单震源的波场的成像结果，而第二项就是多个震源之间震源波场与记录波场间不匹配而产生的交叉干扰像(Dai and Schuster, 2009).若 Γ i和 Γ j存在比较大的差异，即震源i和j距离较远或随机延迟时间差别大，则可以通过叠加和偏移很好地压制交叉干扰像.
由于偏移结果中包含了各个震源波场的贡献，因此从理论上说和常规单震源偏移结果是可以对比的，若不考虑交叉干扰像则直接成像结果和常规单震源成像结果的幅值是一致的.野外实际生产过程中设计观测系统时尽可能把随机延迟时间加大一些、同步激发的震源距离远一些，这样就可以明显降低交叉干扰像，通过偏移成像和叠加又压制了其中部分随机干扰，直接偏移成像可以取得和常规偏移成像方法相差不多的效果.
为了验证上节提出的两种多震源地震数据偏移成像方法的正确性和有效性,本文设计图 4a所示的速度模型进行试验.该模型是中国石化石油物探技术研究院根据实际地质资料建模得到的，模型的浅层中部有一个类似地层尖灭的断层；模型深层是一个背斜构造，利于油气富集与存储.该模型是对实际地质模型的一个简化，可以一定程度上反映断 层和背斜构造在油气勘探中的影响.模型整体呈层状构造，浅层速度小(3000 m/s)，深层速度大(4700 m/s).
模型网格点数：nx=3200、nz=220，网格间距：dx=dz=10 m.试验中的多震源观测系统如图 1所示，单炮数nshot=250，炮间距为80 m，共1600道，道间距为20 m，第一个检波点与第一个震源点的距离xl=6000 m，每一个混合炮的观测范围都等于常规单炮激发总的观测范围.这里取一起激发震源个数N=2，得到125个多震源地震记录，记录长度为8 s，时间采样率为1 ms；震源使用主频为30 Hz的雷克子波.按照上述观测系统进行数值模拟，得到125个含两个震源的多震源地震记录，其中一个如图4b所示.
 | 图4 二维断层-背斜构造模型及其对应的含两个震源激发的混合地震记录(a)速度模型；(b)一个含两个震源激发的混合地震记录. Fig.4 The 2D velocity model with fault-anticline structures and the blended seismic shot gathers with 2 simultaneous sources (a) The velocity model； (b) The blended seismic shot gather with 2 simultaneous sources. |
本文对上述125个含两个震源的多震源地震数据按照文中所述方法进行分离处理，得到了250个单震源地震记录.分离后，各个记录波场都消除了随机延迟激发时间的影响，记录长度减少为6 s.如图4b所示的多震源地震记录分离结果见图 5，可以看到本文方法很好地实现了多震源地震记录的分离，得到了相对应的单震源地震数据，随机延迟激发引起的交叉噪音得到了很好的压制.对分离之后的地震数据进行偏移成像处理，成像结果如图 6a所示，可以看出分离后再偏移成像由于叠前偏移成像过程中的叠加，震源随机延迟激发引起的分离波场中的随机噪音在偏移成像处理中得到了很好的压制，去除了部分随机噪音，整体剖面连续性好、构造清晰.
 | 图5 图4b所示多震源地震记录经过分离得到的两个单震源地震记录 Fig.5 The two separated seismic shot gathers from the blended shot gather that shown in Fig.4b |
为了便于比较，这里对常规单震源地震方法也进行了试验，观测范围和多震源观测范围相同，单炮数nshot=250，震源同样使用主频为30 Hz的雷克子波，记录长度为6 s，时间采样率为1 ms，图6b为常规单震源偏移成像结果，可以看出，多震源分离后再成像结果相比常规单震源偏移成像结果，模型吻合度、构造连续性等方面基本一致，证明本文提出的多震源地震数据分离后再偏移成像方法是正确可行的.
 | 图6 二维断层-背斜构造模型偏移成像结果 (a) 125个含两个震源的多震源地震记录分离得到250个地震记录后再偏移成像结果；(b)常规单震源地震记录偏移成像结果. Fig.6 The migration results of 2D fault-anticline model (a) The migration result of the 250 separated shot gathers from 125 blended shot gathers with 2 simultaneous sources； (b) The conventional migration result by single-source method. |
另外，为了验证激发震源个数对偏移成像的影响，本文在图 1观测系统的基础上，使用5个震源一起激发，得到了50个多震源地震记录，同样完成了对整个测区的观测，其中一个地震记录如图7a所示(从左至右震源的延迟时间分别为：1407、304、1114、1982和1778，单位ms)，可以看出多个单震源记录经过延迟叠加在一起，形成了更为复杂的地震记录.图7b是按照本文方法分离出来的对应第2个震源 的地震记录，与两个震源激发地震记录的分离结果相比，所需要的主要信息已经呈现出来，但有部分残余干扰存在于分离后的记录中；图7c是相应观测条件下只在第二个震源位置激发的常规单炮地震记录，可以看出分离结果还有部分残留，但主要同相轴等信息已经分离出来了.分离后再偏移成像结果如图7d所示，可以看到成像过程的叠加作用压制了分离结果中的部分随机残留干扰，剖面连续性比较好、位置信息比较准确；但也出现了一些大的假象和干扰，主要集中在浅层(图中圈出部分)，这主要是由于分离结果中的残留波场所致；但因为叠前偏移的叠加过程，浅层由于地震波传播距离较小，在成像点上的叠加次数少,而在深层地震波传播距离较大，在成像点上的叠加次数多，这种叠加作用压制了大部分的残留波场.
在主频1.6 GHz的Intel(R) Xeon(R) CPU、内存8G的计算机上进行多震源地震记录分离并进行偏移成像处理，两种观测系统下多震源记录的分离时间大体一样为75 s，偏移成像时间为28773 s，总的计算时间为28848 s.
 | 图7 二维断层-背斜构造模型多震源分离后再偏移成像试验图示 (a)一个含5个震源激发的多震源地震记录；(b)图(a)的分离结果中的第二个震源对应的地震记录；(c) 和图b相对应的常规单震源地震记录；(d)分离得到250个单震源记录的成像结果. Fig.7 The migration result after the blended seismic data separation of 2D fault-anticline model (a) One of the blended shot gathers with 5 simultaneous sources；(b) The 2nd shot gather of the separation results from the data that shown in Fig.a； (c) The conventional single-source shot gather on the location of the 2nd source corresponding to Fig.b; (d) The migration result of the 250 separated shot gathers. |
为了验证本文所述直接偏移成像方法对多震源地震数据的成像效果，这里对数值模拟得到的125个含两个震源的多震源地震数据直接进行偏移成像处理，结果如图8a所示，可以看到这种情况下的直接偏移成像结果比图6a所示的分离后再偏移成像结果信噪比更高，与图6b所示的常规单震源偏移成像精度相当.这说明当同步激发震源个数较少(这 里为2)、激发震源相距较远时，多震源地震直接成 像效果可以取得与单震源地震方法成像相当的精度.
 | 图8 二维断层-背斜构造模型直接偏移成像结果(250个单震源，不同的激发方式) (a) 125个含两个震源的多震源地震记录成像结果；(b) 50个含5个震源的多震源地震记录成像结果. Fig.8 The direct migration result of 2D fault-anticline model (different blended shooting methods with a total number of 250 sources) (a) The direct migration result from 125 shot gathers with 2 simultaneous sources； (b) The direct migration result from 50 shot gathers with 5 simultaneous sources. |
类似地，这里也对如图7a 所示的50个含5震源激发的多震源地震记录进行了直接偏移成像试验，结果如图 8b所示.从成像结果中可以看出直接偏移成像方法对模型的主要构造取得了很好的成像结果，但由于同步激发震源个数的增多以及随机延迟激发时间的影响，整体剖面有一些随机干扰；同样 由于浅层地震波在成像点上的叠加次数少,深层地震波在成像点上的叠加次数多的原因，这种干扰在浅层尤为明显，深层干扰较小、效果较好；整体成像剖面也会出现一些大的假象，在成像剖面中分布比较分散(图中圈出部分).因此，多震源勘探中必须合理地设置观测系统，尽可能使同步激发的震源相隔地远一些，降低相互之间的干扰，提高成像精度.为了更清楚地比较直接偏移成像的交叉噪音干扰，这里取成像结果中2070~2130网格上的成像结果进行对比分析，道集显示如图9所示，(a)为两个震源激发成像结果的道集显示，(b)为5个震源激发成像结果的道集显示，可以清楚地看到随着同步激发震源个数的增多， 成像结果中浅层的干扰也随着增多，深层也有一定的变化，但干扰像明显弱于浅层的.
同样在主频1.6 GHz的Intel(R) Xeon(R)CPU、 内存8G的计算机上进行多震源地震记录分离并进行偏移成像处理，对125个含两个震源的多震源地震数据直接偏移成像处理计算时间为19057 s，对50个含5个震源的多震源地震数据直接偏移成像处理计算时间为7344 s，相比分离后再偏移成像的计算时间28848 s都提高了计算效率，也说明随着震源数目的增加，在提高野外采集效率、降低数据采集成本的同时，室内处理时间明显地减少，降低了处理成本.
 | 图9 图8所示的直接偏移成像结果中第2040—2140网格位置的道显示 (a)两个震源同步激发的情况；(b)5个震源同步激发的情况. Fig.9 The traces from 2040 to 2140 of direct migration result that shown in Fig.8 (a) The condition of blended shooting with 2 sources；(b) The condition of blended shooting with 5 sources. |
Marmousi模型是法国石油研究院推出的模型 数据，是检验波动方程叠前深度偏移成像算法标准 模型数据.该模型构造复杂，可以用来验证本文方法对复杂模型的适应性.
Marmousi速度模型网格大小为nx=243,nz=750，网格间距dx=25,dz=4(单位：m).本试验仍利用如图 1所示的多震源观测系统，取单震源个数nshot=240，炮间距为25 m，接收道数为241道，道间距为25 m，第一个检波点与第一个震源点的距离 xl为25 m，每一个混合炮的观测范围都一样，且等于常规单炮激发总的观测范围.震源使用主频为20 Hz的雷克子波，一起激发的震源个数N=2，激发次数M=120，记录长度为4 s，时间采样率为4 ms.按照上述观测系统进行数值模拟，得到120个含两个震源的多震源地震记录，其中一个(第38个)混合地震记录如图10b所示，按照本文所述方法进行分离处理，得到了240个单震源所对应的地震记录，图10b所示的混合记录对应的分离结果如图10c和图10d所示.可以看到，分离结果中消除了相应的延迟时间，得到了单震源记录的主要信息；由于 波场记录的复杂程度较高，分离结果中还有部分残留.
 | 图10 Marmousi速度模型模拟及其多源数据分离结果 (a)Marmousi速度模型；(b)一个含两个震源的混合地震记录；(c)分离得到的对应左侧 震源的地震记录；(d)分离得到的对应右侧震源的地震记录. Fig.10 Marmousi velocity model and the separation results of the blended seismic data (a) Marmousi velocity model；(b) One blended shot gather with 2 simultaneous sources；(c) The separated shot gather of the left source；(d) The separated shot gather of the right source. |
对通过分离120个混合地震记录得到的240个地震记录进行偏移成像处理，结果如图11a所示.为了便于对比分析，本文还对常规单震源模拟得到的240炮地震数据进行了偏移成像处理，结果如图11b所示.可以看到，分离后再偏移成像的结果构造清晰，和常规地震方法相差不大，本文提出的分离后再进行偏移成像的方法对于复杂模型也是适应的.
应用本文所述直接偏移成像方法，对如图10的 120个含两个震源的混合地震记录进行了直接偏移 成像处理，结果如图12所示.对比图11可以看出,直接偏移成像方法对于复杂模型的成像效果也是不错的，对陡倾等复杂构造都可以很好地成像，本文提出的直接偏移成像方法对复杂模型也是适应的.不足之处是偏移结果中出现了一些随机噪音，信噪比有待进一步提高.
 | 图11 Marmousi模型多震源数据分离后偏移与常规单震源偏移对比试验 (a)分离得到的240炮地震数据偏移结果；(b)常规单震源激发得到的240炮地震数据偏移结果. Fig.11 The comparison of migration after separation and conventional single-source method (a) The migration result from 240 separated shot gathers；(b) The migration result from 240 shot gathers of conventional single-source method. |
为了验证多震源激发成像方法的振幅保真性，这里以图11b和图12所示的偏移结果为例，从以下三点对常规单震源成像方法和直接偏移成像方法的 效果进行对比分析：(1)取偏移成像结果的2450 m位置上的一列数据进行分析.两种方法的幅值变化曲线如图 13a所示，其中蓝线表示常规成像方法成像结果的各个深度的幅值，红线表示直接成像方法成像结果的各个深度的幅值(下同).可以看出两种方法的成像结果数值上有很好的拟合，在一些小的变化剧烈的地方(如图中A、B、C和D点等处)的拟合会有稍微差异.(2)取480 m深度层上x方向各个 位置的数据进行分析.两种方法的幅值变化曲线如图13b所示，可以看出两种方法的成像结果拟合的很好，在振幅突变点(如图中A、B点)也有比较好的效果；直接偏移成像方法在某些地方(如图中C点)会有一定的震荡出现.(3)取1640 m深度层上的各个x方向的数据进行分析.两种方法的幅值变化曲线如图 13c所示，与浅层(上一种情况)相比，两种成像方法的结果整体上拟合的更好，也一 定程度上从侧面说明了深层的成像效果要优于浅层.
 | 图12 Marmousi模型直接偏移成像结果(120个含两个震源的地震记录) Fig.12 The direct migration result of Marmousi model (120 blended shot gathers with 2 simultaneous sources) |
 | 图13 Marmousi模型常规成像方法和直接成像方法结果幅值对比曲线 (a) 2450 m横向位置上的数据幅值对比曲线；(b) 480 m深度位置上的数据幅值对比曲线； (c) 1640 m深度位置上的数据幅值对比曲线.其中蓝线表示常规成像方法，红线表示直接成像方法. Fig.13 The amplitude curves of the conventional migration result and the direct migration result of Marmousi model (a) The curves on the horizontal position of 2450 m；(b) The curves on the depth of 480 m； (c) The curves on the depth of 1640 m. The blue curve represents the conventional migration result and the red one represents the direct migration result. |
通过以上对比分析可知，本文提出的两种对多震源混合地震数据的处理方法对复杂构造模型也是适应的，直接成像方法在一定程度上相比常规单震源地震成像方法是保幅的.
多震源地震勘探方法技术是当今地震勘探的一个热点方向，主要是为了提高勘探效率、降低勘探成本，是一种有着良好发展前景的方法技术，特别是对于三维大面积探区，同步激发震源个数可以设置的更多，道集数也就成倍减少，后续的处理效率也就会有大幅度提升.本文从多震源的本质出发，基于多震源波场的传播规律以及波场的线性叠加原理，提出了多震源数据分离后再偏移成像和多震源地震数据的直接偏移成像方法，两种方法在数值模拟数据的模型试验中均得到了满意的成像效果，表明本文方法是正确的、有效的，可以应用于多震源地震数据的偏移成像处理流程中.主要结论有以下几点：
(1)当多震源中激发震源较少(或相距较远)时，来自各个震源的波场相互干扰相对较弱，分离后再偏移成像和直接偏移成像方法都可以取得和常规单震源地震方法相当的效果，总体精度比较高，这说明本文提出的多震源地震数据的分离后再偏移成像方法与直接偏移方法可以满足叠前深度偏移要求.
(2)当多震源中激发震源较多(或相距较近)时，深层由于叠前深度偏移中叠加次数增多的原因，两种方法偏移成像效果都呈现深层较浅层好的特点；多震源地震记录的分离剖面上残留的一些随机干扰在偏移成像后已经得到了一定程度的压制，但仍有部分假象出现在浅层；直接偏移成像剖面中假象分布比较分散.
(3)多震源地震数据的直接偏移成像方法相比常规单震源地震成像方法是保幅的，特别是在同步震源距离较远时，交叉干扰像会比较弱，成像的保幅性更好.
总的来说，直接偏移成像省去了分离多震源地震数据的步骤，计算效率相比常规偏移成像和分离后再偏移成像有了很大程度的提高；特别是随着多震源中激发震源个数增多，偏移成像的计算效率会更高，但偏移成像效果会有所下降.实际工作中可以把两种成像方法结合起来，直接进行偏移成像处理仍是多震源地震勘探数据处理的一个重要的发展方向；本文是对多震源激发条件的混合波场成像问题进行一些探索研究，多震源地震方法还处于发展阶段，速度建模方法还有待研究，这也是我们下一步的工作方向.
| [1] | 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. |
| [2] | 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. |
| [3] | Bagaini C. 2006. Overview of simultaneous vibroseis acquisition methods. 76th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 70-74. |
| [4] | Beasley C J. 2008. A new look at marine simultaneous sources. The Leading Edge, 27(7): 914-917. |
| [5] | Beasley C J, Chambers R E, Jiang Z J. 1998. A new look at simultaneous sources. 68th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 133-135. |
| [6] | Berkhout A J. 1992. Areal shot record technology. Journal of Seismic Exploration, 1(3): 251-264. |
| [7] | Berkhout A J. 2008. Changing the mindset in seismic data acquisition. The Leading Edge, 27(7): 924-938. |
| [8] | Berkhout A J, Blacquière G, Verschuur D J. 2009a. The concept of double blending: Combining incoherent shooting with incoherent sensing. Geophysics, 74(4): 59-62. |
| [9] | Berkhout A J, Blacquière G G, Verschuur E G. 2008. From simultaneous shooting to blended acquisition. 78th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 2831-2838. |
| [10] | Berkhout A J, Verschuur D J, Blacquière G. 2009b. Seismic imaging with incoherent wavefields. 79th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 2894-2898. |
| [11] | Blacquiere G, Berkhout G, Verschuur E. 2009. Survey design for blended acquisition. 79th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 56-60. |
| [12] | Bouska J. 2010. Distance separated simultaneous sweeping, for fast, clean, vibroseis acquisition. Geophysical Prospecting, 58(1): 123-153. |
| [13] | Chen S C, Cao J Z, Ma Z T. 2001. Prestack depth migration method based on quasi-linear born approximation. Chinese Journal of Geophysics-Chinese Edition (in Chinese), 44(5): 704-710. |
| [14] | Chen S C, Ma Z T. 2006. Generalized synthesis of seismic data and its migration. Chinese Journal of Geophysics-Chinese Edition (in Chinese), 49(4): 1144-1149. |
| [15] | Chen S C, Wang H C, She D P. 2012. Migration of generalized random synthesis of seismic data. Oil Geophysical Prospecting (in Chinese), 47(6): 868-872. |
| [16] | Chen S M, Liu H, Li Y M. 2001. The area shot method for land 3D seismic acquisition. Progress in Geophysics (in Chinese), 16(3): 58-67. |
| [17] | Dai W, Schuster J. 2009. Least-squares migration of simultaneous sources data with a deblurring filter. 79th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 2090-2994. |
| [18] | Hampson G, Stefani J, Herkenhoff F. 2008. Acquisition using simultaneous sources. The Leading Edge, 27(7): 918-923. |
| [19] | Howe D, Foster M, Allen T, et al. 2009. Independent simultaneous sweeping in Libya-full scale implementation and new developments. 79th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 109-111. |
| [20] | 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. |
| [21] | Ikelle L. 2007. Coding And Decoding: Seismic data modeling acquisition and processing. 77th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 66-70. |
| [22] | Howe D, Foser 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. |
| [23] | Mahdad A, Doulgeris P, Blacquiere G. 2011. Separation of blended data by iterative estimation and subtraction of blending interference noise. Geophysics, 76(3): 9-17. |
| [24] | Moerig R, Barr F J, Nyland D L. 2002. Simultaneous shooting using cascaded sweeps. 72th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 74-76. |
| [25] | Rietveld W E A, Berkhout A J, Wapenaar C P A. 1992. Optimum seismic illumination of hydrocarbon reservoirs. Geophysics, 57(10): 1334-1345. |
| [26] | Romero L A, Ghiglia D C, Ober C C, et al. 2000. Phase encoding of shot records in prestack migration. Geophysics, 65(2): 426-436. |
| [27] | Sun P Y, Zhang S L, Liu F. 2002. Prestack migration of areal shot records with phase encoding. 72nd Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 1172-1175. |
| [28] | Ting C O, Zhao W. 2009. A simulated wide azimuth simultaneous shooting experiment. 79th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 76-80. |
| [29] | Wang H C, Chen S C, Zhang B, et al. 2013. Separation method for multi-source blended seismic data. Applied Geophysics,10(3):251-264. |
| [30] | Wu R S, Chen L. 2001. Beamlet migration using Gabor-Daubechies frame propagator. 63rd EAGE Annual Meeting, Extended Abstracts, 1356-1359. |
| [31] | Zhang G, Zhang W S, Hao X J. 1999. Prestack depth migration with common-shot and synthesis-shot records. 69th Ann. Internat Mtg., Soc. Expi. Geophys., Expanded Abstracts, 1469-1472. |
| [32] | 陈树民, 刘洪, 李幼铭. 2001. 适于窄线三维地震资料的面炮方法. 地球物理学进展, 16(3): 58-67. |
| [33] | 陈生昌, 王汉闯, 佘德平. 2012. 地震数据广义随机合成的偏移成像. 石油地球物理勘探, 47(6): 868-872. |
| [34] | 陈生昌, 马在田. 2006. 广义地震数据合成及其偏移成像. 地球物理学报, 49(4): 1144-1149. |
| [35] | 陈生昌, 曹景忠, 马在田. 2001. 基于拟线性Born近似的叠前深度偏移方法. 地球物理学报, 44(5): 704-710. |