干涉测量是以光波干涉原理为基础的一门测量技术，研究成对信号之间干涉现象，获取成对信号之间的差异信息，已广泛用于天文学、测绘学和机械工程等学科，比如光学干涉测量、雷达干涉测量(廖明生和林珲，2003)和全息干涉测量(熊秉衡和李俊昌，2009)等.犹他大学Schuster(2001)首次提出“地震干涉测量”概念(Schuster，2009).地震干涉测量研究地震成对信号之间干涉现象，即两个接收点位置上记录互相关，构建一个新的地震响应，这个响应就像一个接收点上激发、另一个接收点接收的响应.

斯坦福大学Claerbout(1968)导出一维透射响应和反射响应之间关系式，后来命名为声波日光成 像(acoustic daylight imaging).巴黎大学Fink(1997)的时间反转(time-reversal)声学超声波试验，重新 引起众多学科对多次散射尾波的研究兴趣.Schuster(2001)将天然地震中互相关方法延伸到地震勘探领域.Wapenaar等(2002)用格林定理严格证明了Claerbout猜想.Snieder等(2002)用地震干涉测量 法从天然地震尾波中提取了地球信息.Curtis等(2006)指出地震干涉测量法将噪声变成信号.Larose等(2006)和Cao等(2008)研究了多次散射使地震干涉测量具有超级分辨率和超级叠加性能.Snieder等(2008)、Vasconcelos等(2009)和Fleury等(2010)研究了声波散射介质的表定理和格林函数恢复.Ramírez和Weglein(2009)系统介绍了格林定理在地震资料处理中所起的作用.Wapenaar等(2010a，2010b)从物理上由时间反转声学推导了地震干涉测量，解释了地震干涉测量基本概念.

利用地震干涉测量技术，在天然地震中把被动地震测量(如环境地震噪声或微地震响应)变成确定性地震响应，用于反演地壳结构(Bensen et al., 2007)；在地震勘探中主要用水平井中接收数据来消除复杂表层影响(又叫做虚源法)(Bakulin and Calvert， 2006)，或把VSP数据转化为井间数据或反射数据，将多次波转化为一次波，扩大照明区，消除复杂介质影响，改善盐墙、深层构造和复杂构造成像效果(Willis et al., 2006；Vasconcelos et al., 2007；van der Neut et al., 2010；Schuster，2013).在理论研究方面，建立了地震干涉测量与逆散射成像之间联系(Vasconcelos and Snieder， 2008；Curtis and Halliday， 2010；Halliday and Curtis， 2010)；在消除噪声方面，地震干涉测量技术用于压制面波和面波散射波(Dong et al., 2006；Halliday and Curtis， 2008；Halliday and Curtis， 2009；Forghani and Snieder， 2010；Xue，2010；Halliday et al., 2010；Kimman and Trampert， 2010；Duguid et al., 2011；Poletto and Bellezza， 2012).到目前为止，仍未见到关于地震干涉测量技术分离常规山地地震勘探数据中近地表散射波的实例.

地震干涉测量(SI)、逆散射序列(ISS)(Weglein et al., 1997)和表面多次波衰减(SRME)(Verschuur et al., 1992；Berkhout and Verschuur， 1997；Verschuur and Berkhout， 1997)具有如下共同点：(1)都是基于波动理论的克希霍夫积分方法；(2)都不需要地下介质信息；(3)都适用于复杂介质；(4)预测的波场都包含了多个震源子波信息；(5)都属于数据驱动方法，都是用地震数据本身作为算子来预测全部或部分波场；(6)都是由预测和相减两个关键步骤构成，所以它们之间存在必然联系.它们之间差异有：(1)出发点不同，SI是从格林定理出发，综合了波动方程、时间反转和散射理论，ISS是从散射理论出发，SRME是从反馈迭代模型出发；(2)适用性不同，SI适用于陆上和海上地震数据，而ISS和SRME仅适用于海上地震数据，假设自由表面是水面；(3)预测的波场不同，SI可以预测整个波场，包括多次散射波场，而ISS和SRME预测自由表面多次波和层间多次波；(4)数据算子不同，SI用全波场作为算子预测全波场、用直达波场作为算子预测散射直达波、用面波波场作为算子预测散射面波，而ISS和SRME用一次反射波波场作为算子预测多次波；(5)计算方式不同，SI是炮内计算后再求和，ISS是在频率-炮点波数-检波点波数域计算，SRME是在频率-炮点空间-检波点空间域计算，关键区别在于ISS和SRME由于事先不知道一次反射波，而是用全波场代替一次波作为算子，所以采用迭代的方式，当然迭代的同时可以预测高阶多次波，所以SI计算量最小，不涉及大型矩阵计算，不存在空间假频问题，不要求数据规则采样，容易实现三维计算，而ISS和SRME需要共炮集和共接收点道集波场同时参与计算，计算量大，并带来空间假频和规则采样问题，三维计算较难；(6)数学和物理解释不同，用稳相分析可以简单明了解释SI原理，而ISS和SRME难以做到这一点，但可以用SI来直观解释ISS和SRME，例如在预测自由表面多次波过程中ISS和SRME的核心就是褶积型干涉测量，在预测层间多次波过程中ISS和SRME的核心就是褶积型和互相关型干涉测量组合.尽管如此，这三种思路值得互相借鉴.

不均匀复杂近地表严重影响了地震勘探数据的质量，尤其是沿地表传播的直达波、折射波和面波遇到强不均匀体或遇到突变地形时产生的强散射波，这些散射波常常叠置在体波数据中，影响范围宽泛，在时间上、空间上和频率上都与体波叠置，更严重的是散射波能量较强，所以在山地等复杂地表区分离近地表散射波，对提高体波信噪比具有现实意义.本文首先将地震干涉测量基本理论与散射理论结合起来，导出了近地表散射波地震干涉测量表达式，分为互相关型和褶积型表达式，它们由实际波场和背景波场干涉测量构成，为近地表散射波分离提供理论基础.然后用实际山地地震资料说明地震干涉测量法衰减近地表散射波效果，最后给出讨论和结论.

地震干涉测量法代表不同接收点上地震记录互相关产生虚拟震源(简称虚源)地震响应的原理.在任意三维不均匀无损耗介质中，由Rayleigh互易定理和时间反转不变性原理，可以获得地震干涉测量格林函数表达式(Wapenaar and Fokkema， 2006；Wapenaar et al., 2010c)：

(1)式左边G(r_A，r_B，ω)表示待求的频率域虚源格林函数，ω为圆频率，r_B为虚源位置，r_A为接收点位置； R 代表实部，即2 R G(r_A，r_B，ω)=G(r_A，r_B，ω)+G^*(r_A，r_B，ω)，上标*表示复共轭.(1)式右边是一个面积分，积分面S为一个闭合面，所包围的体用V表示， n 表示S的向外单位法向矢量；被积函数中i表示虚部单位，ρ(r)为声波介质密度， r 位于S上，表示震源位置，G(r_A，r，ω)和G(r_B，r，ω)是频率域格林函数，震源位置为 r，接收点位置分别为r_A和r_B，并且r_A和r_B都在V内；为哈密顿算子，表示G的梯度，“·”表示矢量点积.(1)式就是互相关型地震干涉测量的基本表达式，表明两个接收点上格林函数互相关的面积分等于这两个接收点之间因果关系格林函数(Causal Green′s function)和反因果关系格林函数(Anticausal Green′s function)之和. 这里的格林函数不仅包含了直达波，还包含了一次散射波和多次散射波.褶积型地震干涉测量的基本公式为：

为了研究地震散射波表达式，将地震干涉测量基本理论与散射理论结合起来.在散射理论中，将实际声波介质模型分解为背景声波介质模型和扰动声波介质模型，它们所对应的波场分别叫做实际波场、未扰动波场和扰动波场，其中未扰动波场又叫做背景波场(或直达波场)、扰动波场又叫做散射波场(Weglein et al., 2003；Fleury et al., 2010).由此，可 以得到散射波场互相关型地震干涉测量基本表达式：

通常用直达波场代替背景波场，导致直接利用(3)和(4)式计算散射波场的效果差.虽然Ramírez和Weglein(2009)给出的直达波场地震干涉测量表达式(即式(37))，由实际波场与背景波场互相关，可以得到散射波场，但是得到的结果中却多出了一项，那就是背景波场的虚部.而Vasconcelos等(2009)和Fleury等(2010)导出的散射波场表达式中都包含了一个体积分(分别为式(12)和式(31))，这个体积分就是散射理论中常见的Lippmann-Schwinger方程(Weglein et al., 2003)，通常是未知的.即使忽略这个体积分，得到的散射波场表达式为散射波场和背景波场的互相关，而不是实际波场和背景波场的互相关.实际上，上述表达式中的散射波场是指整个散射波场，而不是专指近地表散射波场.

为了研究近地表散射波场格林函数恢复的表达式，将散射理论中背景声波介质定义为均匀近地表和地下构造，扰动声波介质定义为不均匀近地表，如图 1所示.实际声波介质模型声波速度为c(r)、密度为ρ(r)、压缩模量为κ(r)=(ρc²)^－1；背景声波介质模型声波速度为c₀(r)、密度为ρ₀(r)、压缩模量为κ₀(r)=(ρ₀c²₀)^－1，其中 r 表示三维空间位置.

图 1 模型分解示意图将实际声波介质模型(a)分解为背景声波介质模型(b)和扰动声波介质模型(c). Fig. 1 The sketch of model decomposed The actual acoustic media model(a)is decomposed into the background acoustic media model(b) and the perturbed acoustic media(c).

经过附录A中推导，可以得到近地表散射波地震干涉测量表达式，这些表达式面积分中被积函数为实际波场和背景波场的互相关或褶积.用 s表示近地表散射波场格林函数，得到的近地表散射波互相关型地震干涉测量表达式为：

其中利用了远场辐射条件(Wapenaar and Fokkema， 2006)，将实际波场和背景波场的格林函数梯度转化为实际波场和背景波场的格林函数.

用中国西部某地区山地实际地震资料测试地震干涉测量分离近地表散射波能力.该测线是一条宽线，观测系统由三条炮线和三条接收线构成，炮线和接收线重合，测线长度为56 km，横跨平原区、山区和砾石区，高程差达到800 m(如图 2，其中纵坐标为高程，横坐标为米桩号.米桩号间距为1 m，后6位代表测线上桩号，前1位代表测线号).该测线有效炮数为2111炮，每炮有1440道接收，分三条接收线，每条接收线上有480道，中间放炮两边接收，道距和线距都是30 m，炮间距60 m，观测系统参数为 7185-15-30-15-7185，采样间隔2 ms，道长8 s，单线纵向覆盖次数为120次.与Halliday等(2010)试验资料的关键差异之一是他们用单点接收的资料，而这里所用的资料是组合接收的，由36个或72个检波器面积组合.激发震源是炸药震源，在平原区和山 区都是用单井激发，药量为20 kg，平原区激发井深 为16 m，山区激发井深在15.1~71.8 m之间；在砾石区采用组合井激发(1~4口)，药量为3~20 kg，激发井深为5.3~34 m.

图 2 试验测线的高程图 Fig. 2 The topographical map of the tested seismic surveying line

为了适应复杂近地表损耗介质，本文采用反褶积核心算法代替互相关核心算法(反褶积核心算法 参考Wapenaar等(2010b)中式(28)).类似Halliday 等(2010)消除地滚波的用法，图 3显示了褶积和反褶积混合型地震干涉测量的配置.在进行地震干涉测量之前，对原始野外数据进行一些预处理：加载道头、剔除废炮废道、能量均衡、初至波自动拾取等.为了分离直达波场，对全波场用初至波走时，结合时间窗口的方法，保证该时窗内包含直达波和折射波.基本流程是用地震干涉测量分离散射波场，再用多道自适应匹配滤波从原炮集中减去散射波场，得到期望的体波波场(Dong et al., 2006).

图 3 褶积和反褶积混合型地震干涉测量配置图 Fig. 3 Configure for convolution-type and deconvolution-type seismic interferometry

这次测试所用的直达波场时窗长度为380 ms，在图 3所示的地震干涉测量配置图中，褶积和反褶积所用的震源片宽度为300 m，反褶积所用的接收点片宽度为600 m.图 4为地形起伏砾石区炮记录，该炮震源位置在第三条测线上(如图 2中五角星所示)，显示的记录位于第一条测线上，该炮接收点排列位于地形起伏较大的砾石区，尤其380道~480道所在的地表高程明显抬升，所以右侧记录上反向散射波特别明显.该炮记录具有砾石区炮记录典型特征，多组面波及其散射波淹没了体波，面波混响和频散现象非常明显.图 5为图 4相应的地震干涉测量炮记录，右侧线性反向散射波和侧面散射波得到加强，同时还预测了面波及其混响.图 6是用多道匹配滤波法将图 4自适应减去图 5的炮记录，该图表明图 4中的反向散射波和侧向散射波得到明显减弱，面波和面波混响得到压制，突出了炮记录中体波信息.说明这种地震干涉测量法能分离散射波，即使在地形起伏区也有效果.

图 4 地形起伏砾石区炮记录 Fig. 4 CSG in the gravel area with rolling topography

图 5 图 4对应的地震干涉测量炮记录 Fig. 5 The seismic interferometry CSG corresponding Fig. 4

图 6 图 4自适应减去图 5的炮记录 Fig. 6 The CSG is the CSG in Fig. 1 subtracted self-adaptively the CSG in Fig. 5

用同样处理流程和处理参数，对比两种数据体的叠加剖面，一种数据体是没有经历地震干涉测量分离近地表散射波的数据体，其叠加剖面也就是常 规处理的叠加剖面，另一种数据体是经过地震干涉测量分离近地表散射波之后的数据体.图 7a显示了常规处理叠加剖面，图 7b为经过地震干涉测量处理的叠加剖面.对比可以看出，在左部平原区这两个叠加剖面基本相当，这说明经过地震干涉测量处理的资料是可靠的；在中部崎岖山地区，后者效果有所改善(具体分析见局部放大图 8)；在右部砾石区，后者明显改善了叠加剖面的质量(具体分析见局部放大图 9).

图 7 常规处理的叠加剖面(a)与地震干涉测量法消除近地表散射波之后的叠加剖面(b)对比 Fig. 7 Comparing(a)the stack profile using the conventional processing flow with(b)the stack profile using the same conventional processing flow after removing the near-surface scattered wave by seismic interferometry method

图 8显示了图 7山区局部放大，CDP范围为 1400~2000，时间范围为2~6 s，这是典型逆冲构造，地表为崎岖山地.图 8a为常规处理叠加剖面，图 8b为地震干涉测量法消除近地表散射波之后的叠加剖面.仔细对比可以发现，经过地震干涉测量法消除近地表散射波之后，图 8b叠加剖面上同相轴有以下改善：①增强了3~3.8 s、CDP1650~2000之间的背斜构造双曲线同相轴连续性，并且加强了背斜左半支同相轴；②加强了右上角高陡单斜层同相轴的连续性；③对于3~3.7 s、CDP1500~1700之间的复杂交叉同相轴，图 8b加强了高陡同相轴，呈现了完整的蝴蝶结形态的同相轴.

图 8 图 7中山区局部放大图(a)为常规处理叠加剖面，图(b)为地震干涉测量法消除近地表散射波之后的叠加剖面. Fig. 8 Zoomed Fig. 7 on mountain area Figure(a)is the stack profile using the conventional processing flow， and figure(b)is the stack profile using the same conventional processing flow after removing the near-surface scattered wave by seismic interferometry method.

图 9显示了图 7砾石区局部放大，CDP范围为3000~3400，时间范围为2~8 s，这是典型平缓的层状构造，地表为砾石区.图 9a为常规处理叠加剖面，图 9b为地震干涉测量法消除近地表散射波之后的叠加剖面.对比可以发现，经过地震干涉测量法消除近地表散射波之后，图 9b叠加剖面上同相轴连续性得到明显加强，信噪比得到明显提高，具体表现为：①图 9b上2~4 s之间的同相轴清晰可见，而图 9a上只能隐隐约约看到同相轴；②图 9b上4~4.8 s之间出现了同相轴，尤其是4.8 s附近的同相轴连续性好，而图 9a上几乎看不到这些同相轴；③在5~8 s、CDP3100~3300之间，图 9a有一组强的单斜干扰，影响了深层资料的信噪比，而图 9b深层信号信噪比较高，体波信息清晰.

图 9 图 7中砾石区局部放大图(a)为常规处理叠加剖面，图(b)为地震干涉测量法消除近地表散射波之后的叠加剖面. Fig. 9 Zoomed Fig. 9 on gravel area Figure(a)is the stack profile using the conventional processing flow， and figure(b)is the stack profile using the same conventional processing flow after removing the near-surface scattered wave by seismic interferometry method.

总之，由以上实际资料炮集和叠加剖面来看，地震干涉测量法能明显分离近地表散射波，尤其明显改善了砾石区的地震资料品质，提高了体波信噪比. 4 讨论与结论

在山区地震勘探中，起伏地形和不均匀近地表 是地震勘探技术面临的挑战之一.除了复杂地下构 造本身影响外，山地地震资料体波信噪比低的首要原因是近地表影响，这种不均匀近地表产生的散射波影响广泛，分离这种近地表散射波是山地地震资料处理的必要环节.

本文为了分离不均匀近地表散射波，将地震干涉测量理论和散射理论结合起来，导出了近地表散射波地震干涉测量表达式，这些表达式为实际波场和背景波场的干涉测量.经过实际资料试验，表明地震干涉测量法能分离近地表散射波，尤其在砾石区见到了良好效果.该技术的优点在于它能适应于起伏地形和不均匀近地表结构，并且不需要起伏地形和近地表速度信息.该技术不仅能分离测线上散射源产生的散射波，还能分离侧面散射波.

地震干涉测量是一项新兴技术，在地震勘探中应用仍处于起步探索阶段.本文作为地震干涉测量分离近地表散射波的阶段研究成果，仍存在不足，比如复杂山地近地表散射波种类多致使其预测不足、多道匹配滤波自适应减法之后散射噪声消除不足等.文中说明了地震干涉测量理论虽然完备，但实际应用中需要近似，主要近似包括：(1)为了代替波场梯度，采用了远场辐射条件；(2)由于观测系统限制，只能用部分炮点位置代替积分面；(3)由于激发接收条件差异，用实际记录代替格林函数，这要求实际记录波场保持一致；(4)没有考虑震源子波的影响；(5)由于难以确定近地表范围和介质参数，通常忽略公式中介质参数变化的影响.这些近似是影响地震干涉测量预测效果的主要原因，另外复杂山地波场变化大，自适应减法也直接影响了最终应用效果.本文试验表明，分离崎岖山区近地表散射波有待进一步研究，改进方向有以下几点：一是地震干涉测量要求输入一致性的地震数据，需进一步研究复杂山地地震波场一致性校正技术；二是借鉴ISS和SRME中迭代思路来预测多次近地表散射波；三是用震源-接收点干涉测量进一步提高近地表散射波预测效果，四是针对近地表散射波特征进一步研究自适应减法，甚至可以探索直接预测体波而避免减法. 附录A 近地表散射波地震干涉测量表达式推导

Vasconcelos等(2009)和Fleury等(2010)给出了散射波场格林函数表达式的推导，这里仅给出近地表散射波场格林函数褶积型表达式的推导，同理可以得到互相关型表达式的推导.

由扰动介质褶积型互易定理，可以得到扰动介质褶积型表定理.当r_A在体V内时，利用Delta函数积分筛选性质，得

现在关注两种情况，一是r_A和r_B都在体V内(r_A，r_B∈V)，这时(A1)和(A3)式严格成立；二是r_A在体V内、r_B在体V外(r_A∈V，r_BV)，这时(A1)和(A4)式严格成立.这些公式由面积分和体积分组成，其中体积分中涉及到背景介质与实际介质的参数差异.(A2)和(A4)式反映了面积分和体积分之间关系.

对于第一种情况，r_A和r_B都在体V内.在散射理论中，对积分面S上格林函数应用齐次边界条件，即满足下列三种边界条件之一：

(1)索莫菲尔德(Sommerfeld)辐射条件其中

(2)狄利克雷(Dirichlet)边界条件G_0，s(r_A，B，r)=0，;

(3)纽曼(Neumann)边界条件 Δ G_0，s(r_A，B，r)· n =0，， 则(A1)式中面积分为零，得

如果背景介质密度等于实际介质密度(ρ(r)=ρ₀(r))，令散射势υ(r)=ω²ρ(r)(κ(r)－κ₀(r))，则(A5)式简化为

对于第二种情况，当r_A在体V内、r_B在体V外时，将(A1)和(A4)式相加，得

致谢 感谢中国科学院地质与地球物理研究所刘洪研究员的指导和鼓励，感谢胡英博士提供的常规处理，感谢科罗拉多矿业学院CWP提供的SU软件平台.